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ABSTRACT 


Automatic pulse shape control is simulated for the AN/FPN-42 and AN/FPN- 
44A tube type transmitters. A linear, time invariant (LTI) pole-zero model is de¬ 
veloped for each transmitter at a typical operating point using the least squares 
modified Yule-Walker method and Shank’s method. LTI models for a range of op¬ 
erating points are catenated to represent observed nonlinear behavior, and observed 
time variations are added. After these combined models are tested, a linear con¬ 
troller based on the method of steepest descent is implemented. These models, the 
control algorithm and transmitter system details such as power supply droop, dual 
rating and noise are then incorporated into a MATLAB simulation program. 

In a variety of realistic tests the control algorithm successfully shaped the 
Loran-C pulse, except that zero-crossing times were not always in tolerance and the 
algorithm showed a sensitivity to noise. The algorithm controlled Envelope-to-Cycle 
Difference, produced an entire Phase Code Interval of pulses while compensating 
for droop and phase code bounce, and produced a near-optimal transmitter drive 
waveform for the transmitter/antenna system using the dummy load. 
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I. INTRODUCTION 


Modernizing old electronic systems has always presented a challenge to design 
engineers, and the U.S. Coast Guard’s effort to redesign the control system for 
its Loran-C transmitters is no exception. Coast Guard engineers have identified 
commercially*available hardware to replace much of the old control equipment. This 
new equipment will be easier to maintain and operate and will allow more Loran-C 
control functions to be automated. To realize this capability, however, new software 
must be developed to perform each function. This is one of the most challenging 
aspects of the redesign effort. 

One of the most important control functions is shaping the pulse produced by 
older classes of Loran transmitters. A Loran receiver uses the envelope of the pulse 
to identify a standard zero-crossing; if the envelope is distorted, the receiver may 
lock onto the wrong zero-crossing, resulting in a large position error. The software to 
shape the pulse automatically requires a reliable algorithm. In this thesis, a control 
algorithm based on the method of steepest descent is adapted to meet this need. 

In order to test the algorithm fully and to provide a tool for future study, a 
detailed MATLAB computer program is developed to simulate two older transmitter 
classes, the AN/FPN-42 and AN/FPN-44A. With no documentation available on 
the theory behind the design of these transmitters, this is an exercise in system 
identification and modeling. With its wealth of linear algebra and signal processing 
functions, MATLAB is an ideal operating environment for this work. 

Many details of the ’42 and ’44A transmitter systems and of their operation 
affect the shape of the transmitted pulse. To make the simulation a realistic one, 
as many of these details as possible are included. Chapter II gives an overview 
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of Loran-C and provides the background needed to understand these details. It 
explains each of the pulse shape tests found in the Coast Guard’s Specification for 
the Transmitted Loran~C Signal [Ref. 1). In Chapter III, mathematical models for 
the ’42 and ’44A transmitters are developed; in Chapter IV, the control algorithm is 
presented. These models and the control algorithm then form the foundation of the 
simulation program described in Chapter V. Chapter V also includes results from 
a variety of tests performed using the simulation program. Finally, conclusions and 
recommendations for further study are given in Chapter VI. 
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II. AN OVERVIEW OF LORAN-C 

A. LORAN-C IN BRIEF 

1. History 

LORAN, short for LOng RAnge Navigation, is a radionavigation sys¬ 
tem developed during World War II by the famous Radiation Laboratory at the 
Massachusetts Institute of Technology. The first version, called Loran-A, was used 
during the war to guide Allied military ships and aircraft in the North Atlantic and 
Pacific Oceans. By war’s end, Loran coverage extended over most of the areas in the 
North Atlantic and Pacific where U.S. forces operated. Loran-A, with its one to two 
nautical mile (Nm) fix accuracy and its range of 600 to 800 miles, was a significant 
factor in bringing the war quickly to an end and in preventing the loss of aircraft 
because of inaccurate navigation [Ref. 2: p. 153). 

After the war, while Loran-A continued to operate, research began on 
a similar system called the Cycle Matching Tactical Bombing (CYTAC) navigation 
system for the U.S. military. In 1958 the U.S. Coast Guard assumed control of the 
CYTAC system, which was renamed Loran-C. By using a lower frequency band of 
90 to 110 kHz instead of 1.7 to 2.0 MHz as in Loran-A, greater range was possible. 
Also, other technical improvements brought more accurate geographic positioning 
(Ref. 3: p. 2-12]. 

At first, Loran-C was used mainly by the Department of Defense. As 
the number and size of ships passing through coastal U.S. waters increased and as 
several new radionavigation systems were developed, it became apparent that the 
U.S. government should designate one system which it would support. In 1974 the 
Secretary of Transportation adopted Loran-C as the official radionavigation system 
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for coastal U.S. waters with a minimum accuracy requirement of 0.25 Nm and a 
minimum reliability of ninety-five percent of the time in the Coastal Confluence 
Zone (CCZ), essentially the area from the shore out to 50 Nm. By the early 1980s, 
the Coast Guard had phased out the last of its Loran-A stations and had extended 
Loran-C coverage over the entire CCZ [Ref. 4: p. 12]. In 1990, at the request of 
the Federal Aviation Administration, the Coast Guard began a project to extend 
Loran-C coverage from coast to coast in the continental U.S. Today the Coast 
Guard operates Loran-C stations in the U.S.; its territories; and in “host nations” 
such as Italy, Japan and Turkey. In addition, Loran-C stations are operated by other 
nations, such as Saudi Arabia, China and Russia. Figure 2.1 shows the locations of 
Loran-C stations now operating in the U.S. 

2. How Loran-C Works 

Like Loran-A, Loran-C is based on time differences (TDs) between the 
signals of a master station and one or more secondary stations. Beginning with the 
master, each of the stations in the “chain” transmits in turn a sequence of short 
pulses. A receiver located in the chain’s area of coverage measures and displays the 
elapsed time between the signals from each station. The time difference between 
the master and secondary indicates that the receiver is located at some point on 
a hyperbolic line of position. A number of time difference lines of position from 
one baseline (one master and one secondary station) are shown in Fig. 2.2. When 
more stations are added, their time difference lines overlay these and form a grid of 
hyperbolic lines. The secondaries, up to four in number, are designated W, X, Y, 
and Z. Given two or more lines of position, the receiver “fixes” its position at the 
intersection of these lines. 

The chain’s master and the secondary stations repeat the sequence of 
pulses at a fixed rate, according to the Group Repetition Interval (GRI) assigned 
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Figure 2.1: Locations of U.S. Loran stations. 

to the chain when it was first installed. Assigned GRIs vary from 40,000 to 99,990 
microseconds, so the chain’s cycle may repeat anywhere from 10 to 20 times each 
second. By convention in Loran-C, elapsed time is generally described in microsec¬ 
onds or nanoseconds, but not in milliseconds. Each station in a chain transmits 
its own pulse sequence with the same GRI. Progressively longer emission delays, 
with reference to the master, are assigned to each secondary so the signals of each 
secondary arrive in the same order throughout the chain’s area of coverage [Ref. 1: 
p. 2-5]. Emission delays for a chain with three secondaries are diagrammed in Fig. 
2.3. 
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Figure 2.3: Emission delays. 

For consistently accurate emission delays, each station’s repetition rate 
(the time in which its pulse sequence repeats) must be exactly equal to the assigned 
GRI, so that the stations’ pulse sequences do not move relative to each other. To 
ensure this, each station operates three cesium time reference standards (clocks), 
which are constantly compared to each other to check for drift and whose accuracy 
is on the order of 10 -12 seconds. Periodically, the clocks of each station are also 
compared to the master station’s cesium clocks. If the stations’ repetition rates 
are all identical, the control station, with data supplied by two or more monitor 
stations in the chain’s area of coverage, remotely adjusts the emission delay of each 
secondary station’s signal relative to the signal of the master. 
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3. The Accuracy, Reliability and Availability of Loran-C 

Absolute accuracy and repeatable accuracy are two measures of the ac¬ 
curacy of geographic positions obtained by Loran-C. Absolute accuracy is a measure 
of the error between a charted and an observed time difference. Different radio prop¬ 
agation speeds over land and water, inclement weather and other factors change the 
geometry of the grid of hyperbolic lines of position and produce errors. Nevertheless, 
Loran-C meets the minimum accuracy requirement of 0.25 Nm ninety-five percent 
of the time throughout the CCZ. In many areas Loran-C places the receiver within 
0.1 Nm (200 yards) from its true position [Ref. 4: p. 167]. Repeatable accuracy, 
on the other hand, is a measure of Loran-C’s consistency. If a receiver is placed 
at a known position, repeatable accuracy measures the error between two or more 
Loran readings taken at different times. This type of accuracy would be useful when 
returning to a favorite fishing spot or finding one’s home channel entrance in the fog. 
Loran-C’s repeatable accuracy is one of its greatest strengths and is often within 50 
feet [Ref. 5: p. 44]. 

Another strength of Loran-C is its reliability, the percentage of the time 
the master and at least two secondary stations in the chain covering a given area are 
operating correctly. The Coast Guard’s published reliability goal is 99.7%, winch it 
has met consistently [Ref. 6]. 

Signal availability, the percentage of the time a single station operates 
within established tolerances, is the cornerstone of Loran-C’s reliability. The Coast 
Guard’s goal for availability is 99.9%, and it has achieved 99.95% over the years 
[Ref. 7], This corresponds to a little more than four hours per year when the 
average Loran-C station is not providing a reliable radionavigation signal. 


8 








4. Loran-C’s Future 

Loran-C will continue as a vital radionavigation system in the U.S. in 
the near future for several reasons. Loran-C receivers are inexpensive (they start at 
about $450), Loran-C coverage (in the U.S. and in many areas overseas) is extensive 
and reliable, domestic Loran-C users number over one million, and the U.S. federal 
government’s commitment to support it remains firm. According to the Federal 
Radionavigation Plan, the satellite-based Global Positioning System (GPS) and the 
Coast Guard’s Differential GPS program will eventually replace Loran-C, but only 
after several years of reliable operation [Ref. 5: p. 44]. Accordingly, the Coast 
Guard will continue to operate Loran-C in the United States for at least ten to 
twenty more years. 

Currently, the Coast Guard is not involved in any type of Loran other 
than Loran-C. Therefore, throughout the rest of this thesis, general references to 
Loran refer to Loran-C. 

B. THE LORAN-C SIGNAL 

1. The Individual Loran Pulse 
a. General Description 

The Loran pulse is the basic component of the Loran signal. The 
designers of Loran chose to use pulses instead of a continuous wave signal to achieve 
desired range and performance characteristics with less power supplied to the trans¬ 
mitter [Ref. 8: p. 33]. The first 65ps of the Loran pulse, called the leading edge, is 
the only part the Loran receiver uses. This part is specified completely by: [Ref. 1: 

p. 21] 

i(t) = (t — r) 2 e“ 2 *‘“ T ^ 65 sin(0.2;rf + *C P ) r < t < 65 + r (2.1) 

where 
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A is a normalization constant related to the magnitude of the 
peak antenna current in amperes, 
t is time in microseconds, 

t is the Envelope-to-Cycle Difference (ECD) in fts, and 
C p is the phase code parameter: 0 for positive, 1 for negative. 

The first 90ps of the pulse are shown in Fig. 2.4. 


IDEAL LORAN c WAVEFORM 



Figure 2.4: Ideal loran pulse. 

The “tail” of the pulse, also called the trailing edge, is not shown 
in Fig. 2.4. The dynamics of the particular type of Loran transmitter shape this 
part. There are two requirements for the tail of the pulse: it must not generate 
significant frequency components outside the 90 to 110 kHz band, and its amplitude 
after t = 500/is must not exceed a threshold level established for the particular 
transmitter. In other words, one pulse must decay essentially to zero well before 























the beginning of the next pulse in the sequence, and it must be well-behaved as it 
decays. 

The important part of the Loran-C pulse is the third negative-to- 
positive zero-crossing, marked in Fig. 2.4. The receiver uses this “standard'’ zero- 
crossing, also called the 30 microsecond point, to find the elapsed times between the 
master station’s signal and the secondary stations’ signals. The receiver can measure 
these time differences accurately and consistently once it acquires, or locks onto, 
this zero-crossing. In this lock-on process, the receiver first tries to find coherent 
energy at 100 kHz. When it locates a Loran pulse, it measures the amplitudes of 
adjacent pulse peaks. Because the fifth and seventh positive peak amplitudes have 
a unique ratio, the receiver is able to locate the standard zero-crossing which lies 
between them. The receiver sets up a strobed window over the zero-crossing and 
keeps measuring it, thus maintaining “lock” on the signal. If the pulse is distorted 
in some way, the receiver may have trouble maintaining a lock on the pulse and in 
some cases may not be able to lock onto it at all. 

b. Specification: Individual Pulse (Four Tests) 

To minimize the problem caused by distorted pulses, the Coast 
Guard has established a strict specification for the individual transmitted Loran 
pulse [Ref. 1]. This specification defines four measures of Loran-C pulse shape and 
establishes tolerances for them. These four tests compare the measured Envelope- 
to-Cycle Difference (ECD), the half-cycle peak amplitudes (ensemble tolerance), 
the half-cycle peak amplitudes (individual tolerance) and the zero-crossings against 
these established tolerances. This Subsection describes each in detail. 

These four tests use a parameterization of the Loran-C antenna 
current pulse, measured in amperes using a current transformer at the transmitter 
ground return. The parameters consist of the first 13 half-cycle peak amplitudes 
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(normalized so the largest positive value of the pulse equals one) and the first 12 
zero-crossings (in ps, relative to the standard zero-crossing). This parameter choice 
highlights those parts of the pulse most important to the receiver and reflects the 
limitations of signal processing hardware available in the 1950s and 1960s. 

The first three tests apply only to half-cycles one through eight 
where the standard zero-crossing is located. The term “transmitted” pulse refers 
to the current pulse measured at the transmitter ground return, not to the pulse in 
the far field. The terms “assigned” and “ideal” are used interchangeably to indicate 
standard or theoretical values as listed in the signal specification. Similarly, the 
terms “actual” and “measured” are used interchangeably to describe the character¬ 
istics of the real-world Loran signal. 

Test 1: Envelope-to-Cycle Difference (ECD). The Envelope-to-Cycle Dif¬ 
ference is m indication of the position in time of the envelope of the Loran pulse 
relative to the position of the zero-crossings. Figure 2.5 shows the first few half¬ 
cycles of three Loran pulses with ECD values of —5, 0, and +5 ps, respectively. A 
negative ECD indicates that the envelope has been shifted left (or appeared earlier 
in time) relative to the zero-crossings. A positive ECD indicates the opposite. The 
ECD of the Loran pulse may be controlled arbitrarily, within specified limits, at the 
transmitter to obtain a desired pulse shape. 

One problem with ECD is that it changes as the pulse propagates. 
First, when the Loran-C pulse is transmitted, a 90° carrier phase shift occurs by the 
time the pulse has reached the far E-field [Ref. 1: p. 21], resulting in a change in the 
ECD of +2.5 ps. Second, depending on ground or ocean conductivity, ECD continues 
to change as the pulse propagates over the earth’s surface. One model predicts that 
for every 100 Nra the pulse travels over the ocean ECD changes by —0.25/i.- These 
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loran pulses with ecd - -3, 0, +3 



Elapttd Time, microseconds 


Figure 2.5: Loran pulses with three ECDs. 

changes are the result of different propagation speeds of the envelope and of the 
zero-crossings [Ref. 9]. 

A large ECD can wreak havoc on a Loran receiver. It may cause 
the receiver to misidentify the 20ps or the 40/xs zero-crossing as the standard 30ps 
zero-crossing. This results in a large time difference error which could translate 
to a position error of a mile or more. Fortunately, changes in ECD because of 
propagation are predictable with some accuracy, and Loran receivers can be designed 
to compensate [Ref. 9], 

A receiver cannot, however, compensate for unpredictable changes 
in ECD at the point of transmission. Therefore, the ECD of the transmitted pulse 
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is controlled carefully. Each station is assigned a local, or transmitted, ECD value, 
usually zero. The station’s actual transmitted ECD must differ by more than ±0.5 fxs 
from the assigned transmitted ECD. This tolerance is just one sixth of the largest 
ECD difference shown in Fig. 2.5. 

Estimating the ECD of a Loran pulse is a complicated process, but it 
can be done iteratively using the values of the first eight half-cycle peak amplitudes. 
Appendix A outlines this procedure. Once the ECD of the transmitted pulse is 
estimated, an ideal pulse with the same ECD may be generated according to Eq. 
(2.1). The half-cycle peak amplitudes of this ideal pulse are used in the next two 
tests, which apply only for transmitted ECD values of —2.5 to +2.5 fis. 

Test 2: Half-cycle Peak Amplitudes (Ensemble Tolerance). The root- 
mean-square (rms) error between the first eight actual half-cycle peak amplitudes 
and first eight ideal half-cycle peak amplitudes must not be more than one percent 
of the peak amplitude of the puke. Specifically, let S p , p = 1,2, • • • ,8, represent 
the “ensemble” of the first eight half-cycle peak amplitudes of the actual antenna 
current waveform, in amperes, normalized so the largest positive value of the entire 
pulse (usually at, or near, half-cycle number 13) equak one. Let 7 P , p = 1,2, • • •, 8, 
represent the ensemble of the first eight half-cycle peak amplitudes of the ideal 
antenna current waveform, in amperes, normalized in the same way. Then, 

< .01 ( 2 . 2 ) 

Test 3: Half-cycle Peak Amplitudes (Individual Tolerances). In the first 
eight half-cycles of the pulse, the largest difference between the ideal and actual 
half-cycle peak amplitudes must not exceed three percent of the peak amplitude of 






the pulse. In half-cycles 9 through 13, this requirement is relaxed to ten percent: 


14 - s p \ 

< .03 

1 < p < 8 

(2.3) 

14 - s,\ 

< .10 

9 <p< 13 

(2.4) 


Test 4: Zero-crossings. Loran transmitters are extremely narrowband ampli¬ 
fiers designed to resonate at exactly 100.00 kHz. They are usually well tuned 
to this frequency, but instantaneous frequency distortions may exist in the Loran 
pulse, especially in the first two half-cycles. Since a Loran receiver depends heavily 
on the time-domain behavior of the Loran pulse when sampling zero-crossings and 
half-cycle peak amplitudes, any instantaneous frequency distortions in the pulse can 
affect the performance of the receiver. A simple frequency domain spectrum analysis 
of the entire Loran pulse may not adequately detect instantaneous frequency dis¬ 
tortions in the pulse. Therefore, a time domain analysis of instantaneous frequency 
covering the first 13 half-cycles of the pulse is used instead. 

The zero-crossing times and tolerances in Table 2.4 have been estab¬ 
lished for the Loran-C pulse (Ref. 1]. Category 1 tolerances are the most stringent 
and are generally applied to the newer generations of transmitters. Category 2 toler¬ 
ances are more lenient and are usually applied to the older transmitters. Reference 
1 lists exactly which category applies in each test for every station in the Coast 
Guard. 

2. The Loran-C Pulse Group 

a. Format of the Pulse Group 

The Loran signal consists of a group of eight individual pulses trans¬ 
mitted in rapid succession. This increases the average signal power available to the 
receiver [Ref. 8: p. 33]. In addition to these eight pulses, the master station also 
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TABLE 2.1: ZERO-CROSSING TIMES AND TOLERANCES 


Zero¬ 
crossing (ns) 

Time (ps) 

Tolerance (ns) 

Category 1 

Category 2 

5 

-25 

±1000 

±2000 

10 

-20 

100 

1500 

15 

-15 

75 

1000 

20 

-10 

50 

500 

25 

- 5 

50 

250 

30 

standard 

(time reference) 


zero-crossing 



35 

5 

50 

100 

40 

10 

50 

100 

45 

15 

50 

100 

50 

20 

50 

100 

55 

25 

50 

100 

60 

30 

50 

100 


transmits a ninth pulse, which helps the receiver to identify the master. The ninth 
pulse, when “blinked” ON and OFF according to a preset code, allows the master 
station to notify a secondary station in the chain that the secondary is transmitting 
a signal outside of specified tolerances. The Loran pulse decays essentially to zero 
in 500/xs; the pulses in the group are 1000/is apart, except that the master’s ninth 
pulse is transmitted 2000 pa after its eighth pulse, 
b. Multipulse Trigger 

From the receiver’s point of view, the standard zero-crossing pro¬ 
vides the time reference for each pulse in the group and from group to group. From 
the transmitter’s point of view, the time reference is the Multipulse Trigger (MPT). 
When it is time for a station to transmit a pulse group, the Loran timer equip¬ 
ment sends 8 trigger signals, spaced lOOOps apart, to the pulse generator (PGEN) 
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(and one more trigger 2000/zs later in the case of the master station). When the 
PGEN receives a trigger signal, it sends a transmitter drive waveform (TDW) to 
the transmitter and a Loran pulse is produced and is radiated from the antenna. 
When controlling the transmitted Loran signal, then, the MPT is used as the main 
time reference point for the Loran signal. In this thesis, the standard zero crossing 
is used only to perform the zero crossing test on an individual pulse, 
c. Pulse Group Phase Coding 

Another reason for the Loran-C pulse group is to distinguish the 
Loran ground wave from Loran sky waves. As in other low-frequency systems, radio 
waves may take multiple paths to reach a receiver. The groundwave follows the 
surface of the earth while skywaves are refracted and reflected by the ionosphere to 
return to the earth’s surface. Generally the groundwave is used to calculate Loran 
time differences. Therefore a skywave, which has traveled a longer path and is 
thus delayed significantly, represents a spurious signal and may cause a large time 
difference error if interpreted accidentally as the groundwave. To distinguish the 
groundwave from skywaves, pulse group phase coding is used. 

Phase coding is based on the fact that only the skywave undergoes 
a change in phase when traveling from the transmitter to receiver. The ionosphere 
refracts and reflects the pulses in the group and changes their phases by an arbitrary 
amount. The groundwave’s pulse group, on the other hand, arrives at the receiver 
without a phase change (except for the 90° carrier phase shift from the near to 
the far field, which affects all the pulses about equally). Phase coding shifts the 
phases of certain pulses in the group by exactly 180° at transmission, according to 
a standard pattern shown in Table 2.2. 
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TABLE 2.2: LORAN-C PHASE CODES 


Pulse 

Station 

Group 

Master 

Secondary 

A 

++--+-+-+ 

+++++--+ 

B 

+--+++++- 

+-+-++— 


A “+” indicates no phase change, and a ” indicates a 180° phase 
change. Equipped with this expected pattern, the receiver successfully distinguishes 
the groundwave from skywaves. Two successive pulse groups, A and B, are required 
to implement this scheme. This two-GRI transmission sequence, called a Phase 
Code Interval (PCI), repeats constantly. 

d. Transmitter Power Supply Droop 

When e’ ht or nine pulses are transmitted in rapid succession, the 
transmitter’s power supply may not recover fully from pulse to pulse. This problem, 
most prevalent in older transmitters, causes the amplitude of each successive pulse 
in the group to decrease. In general, the first pulse in the group is the largest and the 
last pulse is the smallest. The decrease in amplitude of the smallest pulse relative to 
the largest pulse is called the “droop” and is defined in percent. The Coast Guard 
has established droop tolerances, which are included in the pulse group uniformity 
tests described later in this subsection. 

The practice of dual-rating, explained more fully in Subsection B.4., 
accentuates the droop. A dual-rated station, located between two contiguous Loran 
chains, transmits pulse groups for two chains. Since the power supply of a dual-rated 
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station often has less time to recover than that of a single-rated station, dual-rated 
stations are given larger tolerances for droop. 

e. Specification: Uniformity of Pulses Within a Group 
(Three Tests) 

Of the four tests applied to the individual Loran pulse, and ex¬ 
plained earlier in Subsection B.l.(b), the tests of the half-cycle peak amplitudes 
(ensemble tolerance), the half-cycle peak amplitudes (individual tolerance) and the 
zero-crossings are applied only to pulse one. To measure uniformity of the pulses 
within the group, the test of ECD (as described in Subsection B.l.) is applied to 
each pulse. Two more tests, which examine pulse-to-pulse amplitude differences and 
pulse-to-pulse timing differences, are also applied. 

Test 1: Pulse-to-Pulse ECD Differences. This test reflects in a general way 
how the shape of one pulse differs from the shapes of the others in the group. The 
ECD of any single pulse must not differ from the average ECD of all the pulses in 
the group by more than the tolerances in Table 2.3. 

TABLE 2.3: PULSE-TO-PULSE ECD TOLERANCES 



Category 1 

Category 2 

Single-rate 

Dual-rate 

0.5 ns 

0.7 (is 

1.0/xs 

1.5^s 


Test 2: Pulse-to-Pulse Amplitude Differences. The amplitude of the small¬ 
est pulse in a group must not differ from the amplitude of the largest pulse in that 
group by more than the limits in Table 2.4, calculated as follows: 











TABLE 2.4: PULSE-TO-PULSE AMPLITUDE TOLERANCES, OR 
PERCENT DROOP (D) 



where 

/ p fc max is the value of i(t) at the peak of the largest pulse 
ipjt min is the value of i(t) at the peak of the smallest pulse 

Test 3: Pulse-to-PuIse Timing Differences. Pulses two through eight are 
transmitted at consecutive integer multiples of 1000/iS after pulse one. Relative to 
the standard zero-crossings of pulse one, the standard zero-crossings of pulses two 
through eight must meet the tolerances listed in Table 2.5. 

Pulse nine, which follows pulse 8 by 2000/xs in the master pulse 
group only, is used mainly to identify the master signal and is not used for navigation 
[Ref. 1: p. 2-9]. Thus a tolerance is not assigned to its position in time. 

3. Blink and “Out-of-tolerance” 

Whenever a baseline is not useable for navigation, the first two pulses of 
that secondary station’s pulse group are “blinked” ON and OFF repeatedly (0.25 
seconds on, 3.75 seconds off). The Loran receiver passes along this warning to the 
user, often actually blinking ON and OFF the time difference reading on its display 
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TABLE 2.5: PULSE-TO-PULSE TIMING TOLERANCES 



Category 1 

Category 2 

Single-rate 

Dual-rate 

(p — l)10Q0ps ± 25ns 

(p — l)1000ps ± 50ns 

(p — l)1000ps ± 50ns -f C 

(p — l)1000ps ± 100ns -f- C 


Note: p is the pulse number (2 through 8) of the pulses which follow the first pulse within each 
group. C is 0 for positively phase coded pulses; |C| < 150ns for negatively phase coded 
pulses. The standard zero-crossing of pulse one is the time reference within each group. 

panel. The transmitter station or the control station initiates blink for any of the 
following reasons [Ref. 1: p. 2-8]: 

• Time difference out of tolerance, 

• ECD out of tolerance, 

• Improper phase code or GRI, or 

• Master or secondary station operating at less than one half of specified output 
power, or master station off air (not transmitting a signal at all). 

Automatic alarms at the transmitter station and the control station sound when 
these quantities are out of tolerance. 

In the definition of blink, the four tests of pulse number one and the three 
tests of the entire pulse group explained above are conspicuously absent. There are 
at least two reasons for this. The first is that the control station, with its Loran 
receivers, is monitoring the most important aspects of the Loran signal as far as 
the user is concerned: it ensures that a receiver can maintain lock and that the 
time difference is correct. In this sense, the fine details of the pulse which are the 










subject of these seven tests go beyond the minimum requirements of the Loran 
system to keep a useable baseline. The second reason is that most of the Loran 
control equipment suite was designed and built before modern signal processing 
equipment was available, and consequently these demanding tests are not conducted 
continuously either at the transmitting station or control station. 

Instead, during a one-hour period each day designated for “system sam¬ 
ple,” an operator at each transmitter station manually tests ECD, the half-cycle 
peak amplitudes (ensemble tolerance), and the half-cycle peak amplitudes (individ¬ 
ual tolerance), using an oscilloscope to measure the pulse peaks. He or she then 
enters the values by hand into a computer, which performs the tests and records 
the results. If a failed test is not accompanied by one of the conditions requiring 
blink, station personnel usually do not initiate blink, but instead interpret the test 
as an indication that transmitter maintenance is needed. From time to time, station 
personnel perform all seven tests and several more as well using a portable Loran 
Data Acquisition (LORDAC) unit. They use these results to keep the transmitter 
operating properly, but generally do not initiate blink if a test fails. 

These seven tests thus represent a stricter standard than the conditions 
requiring blink and serve as an early warning of possible transmitter system problems 
which may later require blink. Therefore, a pulse out of tolerance in one of these 
seven tests may still be useable for navigation, but this is not a desired condition. 

4. Dual-rating and Dual-rate Blanking 

As mentioned briefly before, a dual-rated station, located between two 
contiguous chains, transmits pulse groups for two chains. These chains always have 
different GRIs, or rates. Since each chain is independently controlled, dual-rated 
stations are subject to competing, and sometimes conflicting, requirements as the 
pulse groups from the two GRIs periodically overlap in time. Since it is undesirable 
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to transmit part of one pulse group and part of another, the conflict is solved by 
transmitting one and suppressing, or “blanking,” the other. Blanking, which relates 
to the synchronization of two rates, should not be confused with blink, an indication 
of an out-of-tolerance condition. 

Implementing dual-rate blanking is straightforward. A dual-rated tube 
transmitting station’s timing equipment sets up a blanking interval over each pulse 
group, beginning 500/is before the first pulse is triggered and ending 1400/is after 
the last pulse is triggered. The timing equipment tracks the two blanking intervals 
as they move in time. When they overlap, the timer sends MPTs for only one of 
the two rates to the PGEN. 

Two methods are used to decide which rate is blanked when an overlap 
occurs. In priority blanking, the same rate is always blanked, generally the shorter 
one. In alternate blanking, the priority role is passed back and forth between the 
rates at a time interval equal to the length of four times the longer GRI [Ref. 1: p. 
2-9]. 

5. Frequency Spectrum Requirements 

The energy that a station transmits outside the assigned 90 to 110 kHz 
band must not exceed one percent of total radiated energy. Furthermore, neither 
the energy below 90 kHz nor the energy above 110 kHz may exceed 0.5% of total 
radiated energy. 

C. PRODUCING THE SIGNAL 

1. The Loran Transmitter 
a. Types of Transmitters 

As mentioned previously, Loran-C transmitters are extremely nar¬ 
rowband amplifiers designed to resonate at exactly 100.00 kHz. The Coast Guard 





currently operates four types of transmitters, as listed in Table 2.6. The three types 
of transmitters with vacuum-tube power amplifier stages represent three generations 
of tube transmitter technology. The fourth generation, the solid-state transmitter, 
is now the state-of-the-art in Loran-C. 

The solid-state transmitter is superior to the vacuum tube transmit¬ 
ter: it has a cleaner output signal, it has a higher ratio of output power to supplied 
line power, it is more robust, and it requires less maintenance than any other trans¬ 
mitter type. It also has an automatic pulse generating and control system. Many of 
the stations equipped with this transmitter are unmanned and remotely operated. 

For all these reasons, the Coast Guard has considered replacing all of 
its older transmitters with the solid-state transmitter. However, the relatively high 
replacement cost ($2 million to $4 million per station) and the impending closure of 
some tube stations have kept the tube transmitters in operation for the foreseeable 
future. When the last AN/FPN-39 transmitters are removed from service in the next 
year or two, the only tube transmitter classes remaining will be the AN/FPN-42 
and the AN/FPN-44/44A/44B/45. The ‘44 variants and the ‘45 are essentially the 
same transmitter with progressively more power amplifier stages and consequently 
greater output power. The ‘42 and the ‘44A, the subjects of this report, adequately 
represent the remaining tube transmitters, 
b. Transmitter Loads 

Each station has two different transmitter loads: the antenna and 
the resistive dummy load. Several types of antennas are in service, and they vary 
in radiated power and range. The two most common types are the 625-ft and the 
700-ft top-loaded monopoles. The radiating part of these antennas consists of a 
single steel tower and an umbrella-like cap of guy wires leading from the top of 


TABLE 2.6: TYPES OF LORAN TRANSMITTERS 


Transmitter 

Designation 

When 

Designed 

Shape 

Control 

Amplifier 

Type 

Peak 

Power (KW) 

AN/FPN-39 

1950s 

Manual 

tube 

250 

AN/FPN-42 

1950s 

Manual 

tube 

300 

AN/FPN-44 A/45 

1960s 

Manual 

tube 

400/2000 

AN/FPN-64 

1970s 

Auto 

solid-state 

400/800 


the antenna down to anchors arranged on the ground in a circle around the antenna. 
A ground plane consisting of underground copper wires radiating outward from 
the base of the antenna every three degrees forms an electrical mirror image of 
the antenna. The antenna is connected to the transmitter through an impedance¬ 
matching tuning coil. The dummy load, a bank of large resistors, is used to perform 
various tests and maintenance procedures at varying power levels. 

At two of the Coast Guard’s research and training sites an antenna 
simulator is available. Essentially a high-power RLC circuit, the simulator mimics 
the function of a Loran antenna and allows Coast Guard personnel to conduct 
research and testing without interfering with Loran chains operating in the area, 
c. Normal Loran Operating Procedures 

There are two transmitters at each station. One transmitter at a 
time continually radiates a Loran signal using the antenna. This is designated the 
“operate” transmitter. Except during maintenance procedures, the second transmit¬ 
ter is kept in a “standby” status, ready to come on-line should a problem occur in the 
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operate transmitter. Periodically the stanaoy and operate designations are switched, 
allowing technicians to perform maintenance on the formerly operate transmitter. 
The standby transmitter may send pulses into the dummy load at any time with¬ 
out disturbing the operate transmitter and its signals. When transmitter switches 
interrupt Loran-C service for less than one minute, the Coast Guard considers the 
station to be transmitting continuously for availability recording purposes. 

d. Nonlinear and Time-Varying Behavior of Tube 
Transmitters 

This thesis incorporates two important assumptions. First, Loran 
tube transmitters are nonlinear devices, but behave linearly at a given operating 
point. This assumption is examined and supported in detail in the next chapter. 
Second, the transfer functions of the tube transmitters also vary with time. As 
transmitter components, particularly the vacuum tubes in the amplifier sections, 
age over days and weeks, their amplifying characteristics change. When components 
are replaced, small step changes occur to the transmitter’s transfer function. When 
the operate and standby transmitters are switched, the pulse shape control system 
encounters a larger step change in the plant’s transfer function. Loran technicians 
minimize these effects by a great deal of hard work, but the effects still exist to some 
degree. In addition to these internal factors, weather conditions, such as ice forming 
on the antenna and high winds (which distort the shape of the antenna slightly) 
introduce other time variations as well. Thus, from the point of view of a Loran-C 
control system, the transfer function of this plant exhibits both slow changes and 
periodic step changes over hours, days and weeks. In the absence of severe weather 
conditions or component failure, the transmitter may be considered time invariant 
for a period of several hours. This assumption is used also in the next chapter. 
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e. Transmitter Phase Code Balance 


Tube transmitters use a push-pull amplification system, where the 
positive and negative parts of each pulse are amplified by separate banks of tube 
amplifiers. If the transmitter is not balanced properly, the positive half of the signal 
will be amplified more than the negative half, or vice versa. Most often this is 
detected when examining pulses whose phase code is different in GRIs A and B. 
When the pulse flips back and forth, it appears to “bounce.” Phase code balance 
is an adjustment built into the PGEN which increases the magnitude of the TDW 
for negatively phase coded pulses (those pulses which have been inverted by a 180° 
phase change). In this way the phase code “bounce” is removed. 

2. Transmitter Drive Waveforms and Typical Outputs 

A cosine pulse input is used to excite the highly resonant Loran trans¬ 
mitter. A typical TDW and radio frequency (RF) antenna current waveform are 
shown for both the ‘42 and the ‘44A. The terms input and input waveform refer 
to the TDW, and the terms output and output waveform refer to the RF pulse 
captured at the transmitter ground return. Actually both input and output are at 
the same radio frequency. 

In both TDWs, the cosine pulse includes eight full periods or, by Loran 
convention, sixteen half-cycles. To meet spectrum requirements on the ‘44A, a “tail 
drive” circuit adds a damped sinusoid to the end of the input cosine pulse to slow 
the decay of the RF output pulse. This prevents unwanted frequency components 
from appearing in the output. When input half-cycle 16 equals zero, as in Fig. 2.7, 
the tail drive is suppressed. 




3. Controlling the Pulse Shape 


In Figs. 2.6 and 2.7, each input half-cycle has a different peak ampli¬ 
tude. This is the result of the manual control scheme designed for the vacuum tube 
transmitters in the 1950s and 1960s and the pulse generator (PGEN) which imple¬ 
ments it. By turning one of the 16 dials on the face of the pulse generator, the peak 
amplitude of any of the 16 input half-cycles may be adjusted in ten discrete steps. 
The controls of the two PGENs are shown in Fig. 2.8. 

By observing the full-wave rectified RF pulse overlaid with the envelope 
of the ideal pulse, the dials of the PGEN may be adjusted to match the actual RF 
pulse shape to the ideal. The manual control system used for pulse shaping in the 
tube transmitters is diagrammed in Fig. 2.9. 

The manual process of “pulse building” on a tube transmitter is one of 
the most difficult tasks in Loran-C system operation. Adjusting one half-cycle of 
the input affects not just one half-cycle of the output but all of the pulse which 
follows it in time. Also, the discrete steps available on the PGEN may result in 
large jumps in the amplitudes of the output pulse’s half-cycle peaks. Added to this 
are the nonlinearities of the tube transmitters. Even with skilled and experienced 
operators this process can take several hours. Fortunately, time variations in the 
transmitter’s operating characteristics ordinarily change even more slowly, so when 
pulse building, time variations may be ignored. However, because of these slow 
time variations, on each occasion when pulse-building is attempted, the transmitter’s 
operating characteristics are slightly different. From one point of view, this amounts 
to manually controlling in a sixteen-dimensional space a nonlinear device which 
behaves slightly differently each time the control procedure is attempted. 
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Figure 2.6: ( 42 input and output, antenna simulator, pair 30. 
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Figure 2.7: ( 44A input and output, antenna simulator, pair 72. 
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Figure 2.8: Two pulse generators. 



Figure 2.9: Manual control diagram. 
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Understandably, many technicians have opted to control the pulse from 
the “back end” by keeping the transmitter tuned identically at all times instead of 
attempting to compensate for time variations in the transmitter using the PGEN. 
While the results are often more predictable, this approach is certainly time inten¬ 
sive. Although an automatic pulse generation and control system was designed for 
the solid-state transmitter, none has yet been implemented for the tube transmitters. 

D. THE ELECTRONIC EQUIPMENT REPLACEMENT 
PROJECT (EERP) AND THE PURPOSE OF THIS 
RESEARCH 

1. The EERP and its Plan 1 

In answer to the difficulties of manually controlling a tube transmit¬ 
ter and in response to many other considerations beyond the scope of this paper, 
in 1990 the Commandant of the U.S. Coast Guard assigned to the Coast Guard 
Electronics Engineering Center (EECEN) a multi-year project titled the Electronics 
Equipment Replacement Project [Ref. 7). In identifying portions of the Loran-C 
system requiring a redesign effort, EECEN considered the following: 

• The supportability of current and future equipment, 

• The desire to enhance and expand automation, 

• The need to respond to new system requirements, and 

• The desire to remain in close step with technology. 

This process resulted in five major plans. Plan One, titled “EPA/PGEN/- 
LORDAC Redesign,” calls for a redesign of the entire tube transmitter’s monitor 
and control equipment suite, including the Electrical Pulse Analyzer (EPA) and the 
Loran Data Acquisition (LORDAC) equipment. The new control system should be 
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able to monitor and analyze the Loran pulse continuously, or nearly continuously; 
generate and control a Loran-C pulse in tolerance automatically; and record the 
results necessary to build a database of operational history [Ref. 7]. 

2. The VXIbus Based Loran-C Transmitter and Control System 

In 1990 EECEN began to implement Plan One of the EERP by starting 
project W1180, originally titled “EPA/DPA Redesign” and subsequently renamed 
“Timing and Control Equipment (TCE) Redesign.” In a 1991 report titled “The 
VXIbus Based Loran-C Transmitter and Control System,” Taggart and Turban 
describe a prototype control system constructed at EECEN which, it is hoped, will 
perform these functions [Ref. 10]. A simplified diagram of this control system is 
shown in Fig. 2.10. 

The system works as follows: The computer loads an initial TDW into 
an arbitrary function generator (AFG), which sends the TDW to the transmitter at 
each timer trigger. A digital storage oscilloscope (DSO) samples the RF pulse and 
stores it in the computer’s memory. In a closed-loop fashion, the controller then 
computes a new TDW in an attempt to reduce the error between the actual and 
ideal RF pulses. 

Based on the new generation of Automatic Test Equipment known as 
VMEbus Extensions for Instrumentation (VXIbus), this system will be much smaller 
and will have fewer components than the equipment used today. It will give the 
operator new control capabilities over each pulse in the pulse group. This system 
will thus produce a more consistent Loran signal while reducing maintenance. 

Even though the ‘42 transmitter is not included in the EERP (it will be 
phased out in the next few years), it is simpler to operate than the ‘44A and will 
be valuable when developing this VXIbus system. In the end, the VXIbus control 
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Figure 2.10: Diagram of the VXIbus based control system. 

system will operate with the ‘44A but not with the ‘42 [Ref. 7, 10]. For these 
reasons, both the ‘42 and the ‘44 A are included in this research. 

3. Purpose of This Research 

a. Primary goal: A Control Algorithm 

One of the missing pieces in this VXIbus system is a proven al¬ 
gorithm to generate and control the Loran-C pulse shape automatically. Finding, 
developing and testing such an algorithm is the primary goal of this research. This 
paper contributes directly to Phase II of project W1180, Pulse Generator Redesign 
[Ref. 11]. The other phases of Project W1180 are specifically excluded from this 
paper except where they overlap with Phase II. Furthermore, within Phase II, only 
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the pulse-shaping aspects of the PGEN redesign are considered. Precision timing 
between pulses and pulse groups is excluded. 

b. Necessary Tool: A Computer Simulation Program 

Achieving this goal requires a detailed computer program to sim¬ 
ulate the operation of the ‘42 and ‘44A tube transmitters and those parts of the 
VXIbus control system involved in pulse shaping. There are at least three reasons 
for this. First, testing new types of transmitter drive waveforms, especially in closed 
loop control, is safer on a computer than on a 400 kilowatt transmitter. Second, 
working with a computer simulation is much faster, much easier and much more 
convenient both for the researchers in Monterey, California and for the EECEN 
technicians in Wildwood, New Jersey. Third, in this simulation the researcher con¬ 
trols the transmitter completely and can isolate the effects of different transmitter 
and control system factors on an algorithm. In this way control algorithms may 
be tested more thoroughly in simulation than with an actual transmitter. In the 
next chapter, mathematical models for the ‘42 and ‘44A are developed to use in the 
simulation program. 
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III. MODELING THE AN/FPN-42 AND 44A 
LOR AN-C TRANSMITTERS 

A. INTRODUCTION 

Simulating a dynamic physical system requires a dynamic model, a mathe¬ 
matical representation which transforms an input signal into an output signal just as 
the real system does, under the normal range of operating conditions. In this chap¬ 
ter, explicit mathematical models for the AN/FPN-42 and AN/FPN-44A transmit¬ 
ters are developed and tested. First, the modeling approach and data are described. 
Second, the unit sample response of the ’42 transmitter is identified. Third, a lin¬ 
ear, time invariant (LTI) pole-zero model is developed for the ’42; next, observed 
nonlinearities and time variations are added to the model. After the performance 
of this model is tested, the entire process is repeated for the ’44A. 

B. THE MODELING APPROACH 

1. Discrete-Time Representation 

A discrete-time model represents the transmitters in this research, for 
several reasons. First, the VXIbus is a discrete-time system, and discrete-time 
techniques are most easily transferable to it. Second, this field is well-developed in 
the signal processing literature. Third, working with a discrete-time model is most 
convenient. The data is in discrete-time form already and many useful discrete 
signal processing algorithms and computer programs are available. 
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2. Data from the '42 and '44A transmitters 


a. Data Collection 

EECEN provided eighty-six discrete-time input and output data 
sequence pairs for this project, sixty-seven pairs for the '42 and nineteen for the 
'44A. For twenty-seven of these the dummy load was connected to the transmitter 
instead of the antenna or the antenna simulator. Sampled at 10 MHz by a LeCroix 
9410 digital oscilloscope with eight bits of resolution, each sequence is 4096 points in 
length and covers a time period of 409.6 fis. The input sequence, which is the TDW, 
and the output sequence, which is the RF pulse, were sampled simultaneously on 
channels A and B of the oscilloscope. The two signals are synchronized to within 
100 ns, the length of time between adjacent samples. 

The RF pulses were measured at the ground return line to the trans¬ 
mitter using a Pearson current transformer. Both the TDW and the RF data se¬ 
quences were measured in volts, with the input impedance of the oscilloscope set 
to infinity. Accordingly, this simulation uses volts for both TDW and RF pulse. 
Although the RF pulse is customarily measured in amperes, this difference is only 
a scaling factor and does not affect the validity of the simulation. 

These data sequences are essentially the same as those available on 
the VXIbus system, which uses an eight-bit Tektronix digital storage oscilloscope. 
This similarity strengthens the usefulness of this simulation since a control algorithm 
has nearly identical data to work with in this program as well as in the VXIbus 
system. 

b. Effects of Noise and Quantization 

Figures 3.1a and 3.1b show the power spectra of the above two se¬ 
quences. Figure 3.1c is a closeup of this plot. These are the periodogram estimates 
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Figure S.l: Power spectrum of ‘42 with antenna simulator 
TDW and (b) RF pulse. 























31_2=_i . -J-1-1-1---1-I-1 

25 30 35 40 45 50 55 60 


Sompl* number, k 


Figure 3.1c: Closeup of power spectrum, ‘42, pair 30. 

expressed in decibels [Ref. 12, p. 448]. Sample 41 corresponds to 100 kHz and 
sample 2049 corresponds to 5.00 MHz, one half the sample frequency. Figure 3.1c 
also marks the 90-110 kHz frequency band of Loran-C. 

The transmitter’s behavior as a highly narrowband amplifier is ap¬ 
parent in these figures. In this data pair the signal-to-noise ratio (SNR) increases 
from 51 dB to 61 dB from input to output. Here the SNR is defined in decibels as the 
peak signal value minus the average value of the noise found in the highest twenty 
percent of the frequency range of the power spectrum. Underlying this definition are 
the assumptions that the Loran signal power present in this upper frequency band is 
negligible compared to the noise power and that the noise is white or approximately 
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white. Some of the noise that appears in the output is inherent in the transmitter 
system itself, and some is quantization noise added by the sampling oscilloscope. 
Chapter V includes comparison results when each is varied. 

3. Linearity and Time Invariance 

a. Initial Assumption: Linear, Time Invariant (LTI) Within 

Each Pulse 

In this work, a fundamental assumption was made that the tube 
transmitters behave as linear, time invariant (LTI) systems at a given operating 
point, within a limited time interval [Ref. 13, 14]. 

b. Verifying the Assumption 

Testing a system for linearity and time invariance with all possible 
input sequences would be impossible. However, if the assumption is made that a 
system is LTI under certain operating conditions, its behavior may be completely 
represented by a unit sample response sequence, h(n), the discrete version of the 
impulse response. Then the system output can be represented as a linear convolu¬ 
tion, 

!/(") = £ x(m)h(n - m) = x(n) * h{n), (3.1) 

m=—oo 

where x(n) is the input sequence. With its shifting, multiplying and summing op¬ 
erations, linear convolution tests the crucial properties of LTI systems [Ref. 12, pp. 
24-25, 103]. A good match between a “synthetic" output sequence y,(n) (produced 
by convolving actual input x(n) with the proposed h(n)) and ihe actual output 
sequence y(n) strongly implies that the assumption of LTI behavior is valid. From 
this point, an analysis of the system may proceed using mathematical techniques 
applicable only to LTI systems, subject to further validations as more data becomes 
available or the operating point changes. 
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From an analysis of the data pairs available, the operating point is 
assumed to be a function of the shape of the TDW. If this remains unchanged from 
pulse to pulse, the operating point and the unit sample response also stay the same. 
As explained in the last chapter, the limited time interval in which this assumption 
is considered valid is a period of several hours. 

The ’42 model is developed first for a typical operating point before 
examining other operating points. The TDW and RF pulse shown in Fig. 2.6 define 
this typical operating point. The shape of this TDW, including the 80 ps length 
of the significant part, is representative of the TDWs used in operating the ’42 
transmitters in the Coast Guard. 

Analysis revealed that the ’42 transmitter no longer behaves linearly 
when TDWs are applied which excite the transmitter longer than 90 ps. This 
is apparent in the frequency domain. By linear theory, the bandwidth of an RF 
pulse should be no wider than the bandwidth of the input. However, the RF pulse 
bandwidth remained the same when the input signal bandwidth narrowed, and a 
linear model could not be constructed. Therefore, for TDWs in this category, this 
modeling approach is inaccurate and should not be used. 

c. Adapting the Assumption to Nonlinear and Time-Varying 

Behavior 

Implicit in the assumption above is the idea that when the shape 
of the TDW changes, a different unit sample response may be required to represent 
the transmitter accurately. Catenating a number of LTI models to cover a range of 
operating points implements this idea. If the transmitter’s operating points can be 
identified with sufficient precision, and if a unit sample response sequence can be 
identified for each one, then a system which is non-LTI overall may still be treated 
as LTI. 
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4. Time-Domain Pole-Zero Modeling 


a. Single LTI System 

Each unit sample response sequence is by itself a useful mathemati¬ 
cal model of the transmitter. Because each has a finite length, they may be classified 
as finite impulse response (FIR) systems. Modeling an FIR system as an infinite 
impulse response (HR) system is often more efficient, because there are less param¬ 
eters,. It is also more useful, because the roots of the HR parameters have physical 
significance. This is the approach used here. 

An HR system may be represented by a constant-coefficient differ¬ 
ence equation of the form 

y(n) + a\y(n-l) + --- + a P y(n-P) = box(n)-hb!x(n -1) + • • • + fc$x(n -Q ), (3.2) 


or 


a r y = b T x, (3.3) 


where y(n) represents the output sequence and x{n) represents the input sequence. 
Coefficients ai, • • •, ap and &i, • • •, bq are real and constant. Taking the 2 -transform 
of Eq. (3.2) yields the transfer function 


x _ fep + bji ’ 1 + ••• + b Q z- q 
' 1 -f ajz -1 -j-+ a P z~ p 


(3.4) 


Factoring the numerator and denominator of this rational polynomial produces the 
alternate form 

aw = *-/■)(*-/,)■■■(*-/«) ( 3. 5) 

{z - g\){z - g 2 ) ■ ■ • (z - gp) 


where numerator roots /i, • • •, fq are the zeros of the transfer function and denom¬ 
inator roots g\,-",gp are the poles of the transfer function. Complex pole pairs 
represent the natural or resonant frequencies of the system while zeros represent the 
system’s delays, gains, losses and initial conditions [Ref. 15, p. 3]. 
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Pole-zero modeling in the time domain, as used in this thesis, is the 
process of estimating the poles and zeros of an HR system based on the least squares 
criterion. The a coefficients are also known as the Autoregressive (AR) parameters, 
and the b coefficients are also called the Moving Average (MA) parameters. Thus, 
ARMA modeling is another term for pole-zero modeling, which uses both a and b. 

Pole-zero modeling is characteristically applied to random processes 
[Ref. 16]. Most of these sequences represent systems with underlying dynamics of 
a relatively low order, overlaid by random noise. As Section C of this chapter 
describes, the random noise is expressly filtered out to isolate the system dynamics 
of the transmitter. These dynamics, not the random noise, are what the pole-zero 
model is intended to capture. 

b. Catenated LTI Models 

Each model in the catenation may have a slightly different unit sam¬ 
ple response sequence, but each is still an LTI model. Therefore, the catenated model 
may be expressed as a linear difference equation with non-constant coefficients, of 
the form 

y(n) + ai(t, E n )y{n - 1) + • • • + a P {t , E n )y{n - P) 

= bo(t, E n )x(n) + bi(t, E n )x{n - 1) + • • • + b Q (t, E n )x{n - Q) , (3.6) 

where each coefficient is a function of time, t, and of parameter E n , which accounts 
for the nonlinear behavior of the transmitter at different operating points and is 
defined in Section E of this chapter. Both n and t are integer indices of discrete 
time, but they are used differently. Index n, a multiple of the uniform 100 ns 
sampling interval (for N — 1), appears in equations which represent transmitter 
dynamics within a period of 409.6 /is. This period is too short to experience the 
time variations described in the last chapter. Index t, on the other hand, represents 
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the non-uniform time taken by one iteration of the Loran-C pulse shape controller. 
This sampling interval of t is set to four seconds in the simulation program and is 
lengthened slightly whenever the controller skips a blanked GRI. The time variations 
of the transmitter are also indexed along t. Every t d iterations of t, coefficients 
a(f, E n ) and b(t, E n ) are incrementally changed. Using a realistic value of t d = 300, 
the parameters change incrementally every 15 minutes; after a few hours these 
changes may become noticeable. For analysis and testing, the rate of the time 
variations may be increased by lowering t d . In effect this compresses time scale 
t. These time variations are incorporated into the pole-zero model in Section E of 
this chapter. The process of calculating the non-constant coefficients to produce a 
nonlinear, time-varying model of the ’42 transmitter is also described there. 

The following section of this chapter is devoted to estimating the 
unit sample response of the transmitter, h(n), at the typical operating point. Section 
D contains the algorithms which estimate the poles and zeros for this sequence in 
the time domain. 

C. IDENTIFYING THE SYSTEM UNIT SAMPLE 
RESPONSE (’42 WITH ANTENNA SIMULATOR) 

1. Frequency-Domain Deconvolution and its Numerical Problem 
Building on the assumption of LTI behavior within each pulse, an idea 
used previously in pulse-shaping research on the VXIbus system is adapted to es¬ 
timate the unit sample response: frequency-domain linear deconvolution using the 
discrete Fourier transform (DFT). In this operation, the DFT of the output sequence 
is divided, sample by sample, by the DFT of the input sequence. The resulting com¬ 
plex sequence H(k) in the frequency domain may be interpreted as the DFT of the 
time-domain unit sample response h(n). This sequence h(n) is real and is obtained 
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directly from H(k) using the inverse DFT (IDFT). The fast Fourier transform (FFT) 
speeds the computation of the DFT greatly. 

One significant problem exists with this approach: spurious peaks in the 
frequency domain, caused when output DFT samples are divided by input DFT 
samples close to zero. An ideal filter eliminates most of this numerical noise by 
setting equal to zero the elements of H{k) corresponding to frequencies less than 
50 kHz (sample 40) and greater them 150 kHz (sample 123). Because the sequences 
have been zero-padded to twice their original lengths, sample 82 now represents 100 
kHz. However, spurious peaks still exist in the 50-150 kHz frequency band, as Fig. 
3.2 demonstrates. These spurious peaks distort h(n) significantly; unfortunately, 
zeroing frequencies in this range also distorts h(n). A more sophisticated filter is 
required here. 

2. Removing Spurious Peaks with Median Smoothing 

A nonlinear smoothing technique consisting of running medians and a 
lowpass linear filter can remove the spurious peaks in H(k) [Ref. 17]. First, the 
signal is considered to be the sum of rough and smooth parts R[H(k)] and S[H(k)]: 

#(*) = #[#(*)] + £[//(*)]. (3.7) 

The running median Mu[H(k)], simply the median of the U -point sequence H[k — 
M -f 1 ),'■• ,H(k — 1 ),H{k),H(k + 1 ),•••,J/(fc + M — 1), replaces sample H(k). 
Here U is an odd integer and M = {U + l)/2. This smoother separates the rough 
and smooth parts of the signal by removing single samples with large errors. The 
smoother’s output effectively follows a low-order polynomial curve without distort¬ 
ing the surrounding samples as a linear filter would [Ref. 17, p. 158]. This smoother 
preserves sharp discontinuities, a property useful in many applications. 
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Figure 3.2: H(n) for ‘42 with antenna, pair 30, with antenna simulator 
(a) Magnitude and (b) phase. 
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Often, however, a linear filter is used in combination with the median smoother to 
soften these discontinuities. Here a succession of a five-point median smoother, a 
three-point median smoother and a three-point lowpass hanning filter is used. The 
hanning filter has unit sample response 

h(n) = 1/4, n = 0 

1/2, n = 1 

1/4, n = 2. (3.8) 

This combination smoother is applied twice each to the parts of H(k) corresponding 
to frequencies 50-90 kHz (samples 41-74) and 110-150 kHz (samples 91-122). The 
smoother is not used in the Loran frequency band, 90-110 kHz, since no spurious 
peaks were observed in this band. 

This approach proved to be quite effective, as Fig. 3.3 demonstrates. 
The presence of more than one sharp discontinuity within U points can reduce the 
effectiveness of this technique, as samples 114-118 of the phase plots demonstrate. 
The estimated unit sample response, which appears in Fig. 3.4a, represents the 
relatively low-order system dynamics virtually free of random noise and of the nu¬ 
merical errors inherent in the frequency-domain deconvolution technique. This h(n) 
is tested using linear convolution as before, with another similar input and out¬ 
put data sequence pair. Figure 3.4b shows that the estimated sequence correctly 
represents not only the resonances of the transmitter but also the amplitude and 
phase. 

Simulation tests using h(n) from both filtered and unfiltered H(k) show 
that in the absence of spurious peaks in the 50-150 kHz band, the entire filtering op¬ 
eration increases the mean squared error (MSE) between the synthetic pulse and the 
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Figure 3.3: Smoothed H(n) for ‘42 with antenna, pair 30, with antenna 
simulator, (a) Magnitude and (b) phase. 
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ESTIMATED UNIT sample RESPONSE h(n) (PAIR 30) 



(a) 

ACTUAL it SYNTHETIC OUTPUTS ( h(n) WITH PAIR 52 ) 



Figure 3.4: (a) Estimated unit sample response h(n) y pair 30, with antenna 
simulator, and (b) actual and synthetic RF pulses, y(n) and y # (n), LTI 
model, with pair 52. 
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actual pulse by a maximum of two percent. The mean-squared error between two 
arbitrary sequences w x (n) and tt> 2 (n), each of length X, is defined as 

Cm, = 7 £ K (n) - u? 2 (n)] 2 . (3.9) 

L nssO 

When spurious peaks are present in this band, the filtering operation reduces greatly 
the mean squared error and makes an unusable h(n) into a usable one. 

This technique provides a quick and accurate way to estimate the unit 
sample response of the transmitter for any RF pulse, if the TDW is also provided. 
Now a pole-zero model may be constructed for this sequence. 

D. A POLE-ZERO MODEL OF THE SYSTEM UNIT 
SAMPLE RESPONSE (’42 WITH ANTENNA 
SIMULATOR) 

1. Sampling FVequency Considerations 

As mentioned previously, the data sampling frequency /, = 10 MHz is 
quite high relative to the Loran-C frequency band, 90-110 kHz. Ideally, a lowpass 
filter with cutoff frequency f c = 110 kHz could be applied and the data could be 
sampled at /« = 220 kHz without losing any significant Loran-C information. Thus, 
from one point of view, the data has been oversampled by a factor of 45. 

EECEN personnel sampled the data at /, = 10 MHz to provide the most 
information possible for this research. In particular, the high f, selected allows a 
more thorough analysis of the system noise and provides accurate zero crossing 
times. The push-pull amplification of the tube transmitters may cause zero-crossing 
distortion from time to time, so this extra information is valuable. 

If desired, a lowpass filter may be applied to these data vectors and 
they may be resampled at a lower rate (i.e., decimated) for analysis and simulation. 
In fact, many advantages exist in this approach: the data vectors are shorter and 
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require less storage; the speed of the modeling and simulation programs increases; 
the poles and zeros are not as close to the real axis and to the unit circle in the 
z-plane, yielding a more stable system; and the modeling algorithm performs better 
when the frequencies of the roots are farther apart from each other. 

Disadvantages also exist in decimating these vectors, however. In the 
presence of quantization and other noise, a great deal of resolution in the zero¬ 
crossing times is lost. For example, at /, = 1.25 MHz (corresponding to a decimation 
factor of N = 8), the maximum error allowed for the ’44A pulse’s 40 ps zero¬ 
crossing is 50 ns, one-sixteenth the sampling interval. Zero-crossing times estimated 
by interpolation at this f t are not as accurate as when interpolated at /, = 10 
MHz. Also, for sampling frequencies less than 10 MHz, interpolation is necessary 
when estimating the half-cycle peak amplitudes. This is because the samples do not 
fall at exactly the peak of each half cycle in general. This interpolation introduces 
noise which may cause problems in closed-loop control. At f t — 10 MHz the peak 
estimation error is less than 0.1 percent of the peak value and may be safely ignored. 

To reflect these two competing criteria, the data was analyzed at four 
different sampling frequencies: 1.25 MHz, 2.5 MHz, 5 MHz, and 10 MHz, corre¬ 
sponding to decimation factors N = 8, 4, 2, and 1, respectively. The best overall 
performance occurred at 10 MHz, and so the following sections on pole-zero model¬ 
ing are presented at this sampling frequency. 

2. Technique for Estimating the AR Parameters: The Least Squares 

Modified Yule-Walker Method 

A number of techniques for linear modeling are based on the statistical 
characteristics of the signal being modeled. In this section, the least squares modified 
Yule-Walker method is used to find the a parameters of the HR model of h(n). 
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The autocorrelation function of h(n) is defined as 

OO 

Rh(i) = 51 h(m)h(i + m), -oo<i<oo. (3.10) 

m=-oo 

From Eq. (3.2), Rk{i) can be expressed in the difference equation form 
Rh(*) + &iRh{i — 1 ) +- 1 - apRh(i — P ) 


= boh(i) + bih(i — !) + ••• + 6qh(t — Q ), (3.11) 


which can be written in matrix form [Ref. 16, p. 565]: 


r 1 


’ 7 ' 
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Here R fl has dimensions (Q -f 1) x (P + 1) 
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with L > P + Q. The components of vector 7 are given by 


(3.12) 


(3.13) 


(3.14) 


70) = 5Z b(j)h m (m-j). 


(3.15) 


with b(j) defined as 

^ = { 0; otherwise^ ' < 316 > 

The lower partition of Eq. (3.12) is solved first to yield an estimate of a. If the 
correlation function and the model order P , Q were known exactly, only P equations 
would be required to find a, and Rg would need only P rows. The remaining 
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L — (P + Q) rows of Re would be redundant. However, because these quantities 
are not known exactly, the overdetermined set of equations is more appropriately 
solved for a in the least squares sense. Let t be the error vector that results from 
an arbitrary choice of a: 

Re* = e . (3.17) 

The solution of the following equation minimizes e: 

(RfRc)* = [ r . (3.18) 

This equation is solved by partitioning Re as 



(3.19) 


and estimating a using the pseudoinverse: 

a = —R£r 0 • (3.20) 

The MATLAB left division command (“\”) provides a method for computing the 
pseudoinverse of a rectangular matrix with a high degree of numerical precision. 
This algorithm is based on the QR decomposition [Ref. 18]. The Singular Value 
Decomposition (SVD) could not be used here because of the large size of R £ (with 
4093 rows in Re, the SVD unitary matrix U is (4093 x 4093) and requires 134 MB 
of computer memory). Results obtained with the SVD using smaller portions of 
Re proved to be less accurate than those obtained with the MATLAB left division 
command when using all of Re- 

3. Technique for Estimating the MA Parameters: Shank’s Method 
If the above statistical approach was continued, vector b could now be 
solved by first calculating 7 , using 

7 = R B a (3.21) 
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and then applying spectral factorization techniques. However, a better time-domain 
match is obtained using the deterministic approach of Shank’s method [Ref. 17, pp. 
510-512, 558-560]. 

Shank’s method begins with the estimate of a found by one of the least 
squares methods, as in the previous Subsection. This all-pole model may be ex¬ 
pressed by the transfer function 

. (3-22) 

where A(z) is the denominator of Eq. (3.4). The desired HR model transfer function 
is then 

H(z) = B(z)H a (z). (3.23) 

Using the all-pole model’s unit sample response M n )> which is derived from H A (z), 
the time-domain modeling error of the pole-zero model is 

e B (n) = h(n)-h A (n)*b(n). (3.24) 

Figure 3.5 is a schematic representation of Eq. (2.4). B{z) is chosen so that the 

sum of squared errors is minimized: 

S«=EM»)f. 0-25) 

n=0 

Then vector b satisfies 

H^b = h (3.26) 

in the least squares sense, where 

MO) o ••• o 

Mi) MO) ••• o 

M<?) M<?-i) ••• MO) 

h A (L — 1) h A (L- 2) k A (L-Q- 1) 
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Figure 3.5: Diagram of Shank’s method. 


and 

MO) 

Mi) 

h = HQ) ‘ 
.HL- l). 

Vector b is estimated using the pseudoinverse, as before: 


(3.28) 


b = HJh. 


(3.29) 


and 


4. The Pole-Zero Model 

By trial and error, model order P = 4, Q = 3 was chosen. Vectors 


a = 


b = 


1.0000 

-3.9856 

5.9645 

-3.9723 

0.9934 

0.0513 

-0.1508 

0.1640 

-0.5650 


(3.30) 


(3.31) 


model h(n) of Fig. 3.4a with the minimum mean squared error. Here a has the form 


[1,01,0], • • • ,ap ]' and b has the form [&o> * • , &q]'. The process of selecting the 

model order is examined in detail later in this section. 
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The poles and zeros of this model are calculated from a and b: 


' 0.9983 e +j00675 ' 
0.9983 e~ i0067i 
poles ~ 0.9984 e +J° 0569 
. 0.9984 e~i° 0569 . 


‘ 1.0361 e+> 00677 
zeros = 1.0361 e~*° 0677 

_ 1.0260 e j0 


(3.32) 


When the elements of a and b of are substituted into Eq. (3.2), a unit sample input 
yields the model sequence fuB(n). Figure 3.6a contains a z-plane plot of these poles 
and zeros while Fig. 3.6b is a time-domain plot of h(n ) and /uafa)- 

Overall, the time-domain match is excellent, indicating that the pole- 
zero modeling algorithm has performed well. This is a non-minimum phase system 
and therefore cannot be inverted because that would result in an unstable system. 
Controlling this system using certain algorithms is now potentially more difficult. 

5. Two Criteria for Selecting Model Order 

The competing criteria of accuracy and simplicity are used to select the 
HR model order. The criterion of accuracy is expressed by two time-domain mea¬ 
surements. The first is the mean-squared error between h(n) and h^B (n)- The 
second is the mean squared error between actual and synthetic RF pulses y(n) and 
y,(n) = x(n) * h(n), where x(n) is the actual TDW sequence corresponding to y(n). 
This is the same simulation test described previously. In this case, however, both 
y(n) and y t (n) are normalized so that the maximum positive amplitude of each 
equals one. This quantity, called the normalized mean squared error (NMSE), mea¬ 
sures the effectiveness of the modeling algorithm by comparing shape and phase 
information while ignoring any difference in the maximum pulse peak amplitudes of 
y(n ) and y,(n). The reason for ignoring the amplitude difference lies in the data. 
The overall transmitter gain for data pairs obtained weeks and months apart was not 
generally the same, perhaps because of the components periodically replaced over 
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Amplitude (volts) 


POLE/ZERO PLOT (PAIR 30) 



(a) 



0 500 1000 1500 2000 2500 3000 3500 4000 4500 

Somple number 

(b) 

Figure 3.6: (a) Pole-zero plot of ‘42 LTI model, pair 30, with antenna 
simulator, and (b) original and model sequences, h(n ) and h A B(n). 
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that time period. Therefore, differences in the amplitudes of the h(n) sequences for 
these pairs may be excused. Special data pairs were obtained to map the relation¬ 
ship between input and output maximum positive amplitudes, and the simulation 
program uses these to scale the output. Thus this problem is not a serious one. 

The criterion of simplicity indicates that a lower model order is better. In 
the simulation program, assigning more poles and zeros takes more time. Therefore, 
increasing the model order without obtaining a corresponding increase in accuracy 
is undesirable. Also, when the model order is unnecessarily high two negative ef¬ 
fects may occur. The first is that the poles and zeros may not be consistent from 
one h(n) sequence to the next. For example, one h(n) may have a complex zero 
pair and a real zero, while the next may have three real zeros. This hampers the 
implementation of the nonlinear model described in the next section. The second 
is that the effective rank of R.£ may be less than P, or the effective rank of H a 
may be less than Q. This may cause numerical problems in the modeling algorithm 
when computing the pseudoinverse. Other indications that the order is too high are 
pole-zero cancellations (when poles and zeros migrate to the same locations and, in 
the transfer function, cancel each other out) and large negative real zeros. 

Selecting model order, then, is necessarily a somewhat subjective process. 
Table 3.1 lists the two criteria and associated remarks for a range of model orders for 
the ’42 transmitter. Pair 30 provides the sequence h(n), as before; pair 52 provides 
the test TDW and RF pulse. Orders below P = 4, Q — 1 were wholly inadequate. 
Model order P — 4, Q = 3 was chosen according to these criteria. AR models 
obtained by the least squares modified Yule-Walker method are not effective at this 
sampling frequency, but at lower sampling frequencies their accuracy approaches 
that of the ARMA models. However, they still require nearly double the number of 
parameters. Perhaps a more deterministic AR modeling algorithm such as Prony’s 


59 










TABLE 3.1: BASIS FOR SELECTING MODEL ORDER, AN/FPN-42 
TRANSMITTER 


p 

Q 

Measure 1 
MSE 

Measure 2 
NMSE 

Remarks 

4 

l 

2.1069 x 10-“ 

1.2325 x 10“ 3 


4 

2 

7.6821 x 10" 5 

1.2633 x lO' 3 


4 

3 

1.8948 x 10-* 

9.0733 x lO -4 

Best overall ******* 

4 

4 

1.8314 x 10- 5 

9.1274 x lO" 4 

4th zero: at z = —350 + jO 

5 

3 

3.2620 x lO' 5 

8.8762 x 10" 4 


5 

4 

1.8585 x 10' 5 

8.9258 x 10' 4 


5 

5 

1.9309 x 10~ 5 

7.4982 x 10‘ 4 

Mtx close to singular 

6 

3 

5.5609 x 10" 5 

9.7328 x lO -4 


6 

4 

1.8415 x 10“ 5 

9.1727 x 10" 4 


6 

5 

4.6451 x 10- 5 

7.3653 x 10" 4 

Mtx close to singular 

6 

6 

1.1680 x 10" 4 

1.3550 x 10“ 3 


10 

0 

2.9922 x lO" 3 

9.7445 x 10' 2 

AR models 

18 

0 

1.2797 x lO' 3 

4.5499 x 10” 2 


24 

0 

2.3607 x lO' 3 

7.3931 x 10~ 2 



Method would produce better AR models [Ref. 15, pp. 88-89; Ref. 16, p. 550]. 
However, that is not the subject of this thesis. Completely deterministic ARMA 
modeling (for example, using Prony’s Method to find a and Shank’s method to find 
b) is not quite as effective here as the statistical/deterministic combination of the 
least-squares modified Yule-Walker method and Shank’s method. 
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E. NONLINEAR, TIME-VARYING MODEL OF THE 
AN/FPN-42 TRANSMITTER 

1. Representing Nonlinearities by Moving Poles and Zeros 

a. Changes in the Positions of Poles and Zeros Caused by 

Changes in TDW Shape 

The transmitter’s unit sample response changes slightly as the shape 
of the TDW changes. The pole-zero models of these sequences are correspondingly 
different also. The pole-zero scatter plot of five data sequence pairs in Fig. 3.7 
illustrates this. All five were obtained within a period of three hours, avoiding time 
variations in the transmitter. The length of time each TDW excited the transmitter 
ranged from 5 pis to 80 pis, which provides a range of differently shaped TDWs. 
The average MSE between h(n) and hxB{n) for these five pairs is 3.4641 x 10 -6 , 
indicating an excellent match. This validates the assumption of LTI behavior at 
operating points other than the typical one described previously. 

b. Assigning Poles and Zeros by Parameter E n 

The apparent trajectories of the poles and zeros in Fig. 3.7 imply 
that the transmitter may be simulated effectively by assigning the poles and zeros 
of the model based on the shape of the TDW. In forming this catenated model a 
reliable way is needed to relate the changes in TDW shape to the trajectories of 
each pole and zero. 

The energy of the normalized TDW (with the TDW’s maximum 
positive amplitude equal to one) can be used to assign poles and zeros according to 
TDW shape. This energy, in units of watt-seconds, is defined as 


x(n) T x(n)T 
{max [x(n)]} 2 fl ’ 


(3.33) 
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Figure 3.7: Pole-zero scatter plot for *42 with antenna simulator for five 
LTI models. 

where sampling interval T is in seconds and load resistance R is normalized to 
one. The choice of parameter E n reflects an assumption that the most important 
characteristic of the TDW shape is the length of time the TDW provides a significant 
level of excitation to the transmitter. As observed before, transmitter behavior 
changes significantly for TDWs that excite the transmitter longer than 100 p.s, and 
the transmitter no longer behaves linearly for TDWs that excite the transmitter 
longer than 90 /zs. It seems reasonable, then, to assume that the length of time 
the TDW excites the transmitter is an important indication of transmitter behavior 
also when this length is less than 90 /zs. If the pole and zero locations are consistent 
from one h(n) to another, if the poles and zeros form a fairly predictable trajectory 
as they move with respect to E n , and if the synthetic RF pulse y,(n) matches the 
actual RF pulse y(n), the approach using E n may be considered valid. 
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Poles and zeros are assigned using E n as follows. First, the trajec¬ 
tory of each pole and zero is assumed to be smooth and continuous, rather than a 
series of discrete steps. The magnitude and phase of the kth root from each data 
pair are plotted, so that there as many data points on the magnitude plot and the 
phase plot of the ibth root as the number of data pairs. A polynomial curve is fitted 
to the points on both plots, and the process is repeated for each root. 

A special MATLAB function M0D2 was written to fit these curves. 
The user selects the order of the polynomial curve and may even add points at his 
or her discretion to stabilize or bend the curve. The original points are denoted “o” 
and any added points are denoted When the user is satisfied with the curve, he 
or she selects the minimum and maximum values of E n and so defines the range over 
which the curve is valid. The program stores these endpoints and the polynomial 
values at the endpoints. These are used when assigning poles and zeros if E n falls 
outside the endpoints. The coefficients Ct, c/_j, • • •, c t ,cq of the polynomial 

ct{E n Y + H-h c\E n + co (3.34) 

are stored as well. The curves fit by this program for the magnitude and phase of 
the inner pole appear in Fig. 3.8. If the pole and zero locations in one h(n) sequence 
are not consistent with those of other h(n) sequences, then that h(n) may create a 
spike in one or more of the parameter curves. These outlier h(n) sequences may be 
discarded to fit the curve more accurately. 

When a TDW is presented to the MATLAB function XMTR, which 
simulates the transmitter, the function first computes E n for the TDW. The func¬ 
tion recalls in turn the stored coefficients of each polynomial, evaluates them at this 
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Figure 3.8: Fitted curves for inner pole pair, ‘42 with antenna simulator 
(5 pairs), (a) Magnitude, and (b) phase. 


































E n , and combines these magnitudes and phases into poles and zeros. Vectors a and 
b are then formed from these roots. Next, a unit sample and Eq. (3.2) produce 
A>is(«), which is then convolved with x(n) to yield y # (n). Finally, the amplitude 
of y,(n) is scaled appropriately, according to the input/output energy relationship 
from special data pairs obtained for this purpose. 

c. Performance of the Catenated Model 

The catenated model, as implemented in the function XMTR, was 
tested using the five data pairs described previously. As an example, the normalized 
actual and synthetic RF pulses y(n) and y,(n) for pair 14 appears in Fig. 3.9. The 
sequences are normalized because of the overall transmitter gain difference discussed 
previously. The match is an excellent one; the plots of the other four pairs of y(n) 
and y,(n) show about the same level of performance. The NMSE values between 
y(n) and y,(n) for all five pairs appear in Table 3.2. As a basis for comparison, 
the same tests are performed using the linear model of pair 30 developed earlier 
in this chapter instead of the catenated nonlinear model. Table 3.2 also includes 
these results. The simulation error of the nonlinear model is about an order of 
magnitude smaller in these tests than the simulation error of a simple linear model, 
for an average reduction of 93.4 percent. For data pair 14 the performance of both 
models is about the same. These results validate this method of representing the 
transmitter’s nonlinearities demonstrate that the model accurately reflects the 
transmitter’s behavior over a wide range of operating points. This extra degree of 
accuracy could make the difference between developing an algorithm that works on 
the nonlinear transmitter and VXIbus control system and developing an algorithm 
that does not. Next, the time variations of the transmitter are added to this model. 
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1 


actual & synthetic PULSES (PAIR 14) 



Figure 3.9: Actual and synthetic RF pulses, y(n) and catenated 

model of ‘42 with antenna simulator. 


Table 3.2: PERFORMANCE IMPROVEMENT IN SIMULATION FOR 
NONLINEAR MODEL VERSUS LINEAR MODEL 


Data 

Pair 

NMSE 

(Linear model) 

NMSE 

(Nonlinear model) 

22 

1.2698 x 10- J 

5.7835 x lO' 4 

23 

3.2698 x 10- 3 

2.3224 x 10” 4 

12 

4.6052 x 10- 3 

1.3439 x lO" 4 

13 

5.6114 x 10~ 4 

8.6276 x lO" 5 

14 

3.9731 x 10‘ 4 

3.9176 x 10~ 4 

Mean 

4.3062 x 10’ 3 

2.8460 x 10‘ 4 

Values 
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2. Representing Time Variations by Changing Polynomial 

Coefficient Co 

a. Changes in the Positions of Poles and Zeros Caused by 

Time Variations 

Pole-zero models were produced for three data sequence pairs with 
nearly identical TDWs, obtained over a period of nine months. Figure 3.10a shows 
the first 1000 samples of the TDWs and Fig. 3.10b displays the pole-zero scatter 
plot of the three pole-zero models. From Fig. 2.10b, it is apparent that time 
variations in the transmitter’s transfer function manifest themselves in the same way 
as nonlinearities: by movements in the pole and zero locations of the HR model. 
These variations are about one-fourth the size of the variations due to changes in 
TDW shape. From what is known about the transmitter’s time variations, the 
poles and zeros can be expected to drift slowly over hours, days and weeks. When 
a different transmitter is switched on line, the pulse shape controller sees a step 
change in pole-zero locations. 

b. Moving Poles and Zeros Using Coefficient co 

Changing coefficient Co slightly in each polynomial as a function of 
time t (every tj iterations) is an effective way to model slow time variations in the 
transmitter. As a result, each polynomial curve drifts up and down independently 
of the others with respect to time t. This causes the poles and zeros to drift also 
with respect to t. Introducing a larger random change in the Co of each polynomial 
implements a transmitter switch; the curves, and the poles and zeros, jump to a new 
location. The size of the change to Co is different for each curve, characteristically 
small for the poles and larger for the zeros. Initially the maximum allowable change 
of Co can be set to approximately one fourth the difference between the maximum 
and minimum values of the polynomial curve in the valid range of E n . 
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Sample number, n 



Figure 3.10: (a) Three TDWs (from pairs obtained over a 9-month period) 
for *42 with antenna simulator, and (b) pole-zero scatter plot. 
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c. Simulating Slow Drifts Using an AR2 Random Process 
The following AR2 random process driven by white noise produces 
slow drifts generally in the range (—1,1): 

d(t) - 2(0.996)d(t - 1) -I- (0.996) 2 d(t - 2) = w{t ). (3.35) 

Here d(t) is the output and w(t) is a white noise input with variance *.25 x 10~ 8 . 
The system’s double pole at 0.996 filters the noise and produces the slow drifts. The 
exact value of these parameters were chosen by trial and error. An example output 
d(t) over 1000 iterations appears in Fig. 3.11. Chapter V contains details on how 
parameter drift is implemented in the simu- ition program. 


SLOW 0RIFT d(t) 



Figure 3.11: Slow drift produced by AR2 process. 

3. The Combined Nonlinear, Time Varying Model 

The catenated model from Subsection E. 1. and the time variations 
described in E. 2. combine to form an accurate nonlinear, time-varying model of 
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the ’42 transmitter. Only one more modification is necessary: the inclusion of all 
the data pairs available. Figure 3.12 shows the magnitude and phase of the inner 
pole (positive phase) using nineteen data pairs for the ’42 transmitter. Appendix C 
contains a complete set of these curves. As mentioned before, on occasion the user 
may elect to remove one or more outliers; therefore all the points in the five-point 
curves of Fig. 3.8 may not necessarily appear in these final curves. 

From these curves, new maximum values of Co may be determined sim¬ 
ply by estimating a standard deviation of the data points by eye from the curve. 
Granted, some of the deviation may be from nonlinearities as well as time vari¬ 
ations. If, however, the curve is allowed to drift to cover most of the points (it 
will drift approximately one standard deviation from its normal position), the com¬ 
bined model effectively duplicates nearly all the behavior observed in the data, both 
nonlinearities and time variations. 

The combined model thus meets the need for which it was designed: it 
simulates all the observed nonlinear and time-varying behaviors of the transmitter 
during the convergence of a pulse-shape control algorithm. As the control algorithm 
changes the TDW shape, the model’s transfer function changes also, just as the real 
transmitter’s transfer function does. This dynamic feature is essential for testing 
control algorithms. One of its disadvantages, however, is that the simulation errors 
of the combined model are now higher than those of the linear model. Tested with 
all nineteen data pairs, the linear model’s average NMSE is 2.8531 x 10 -3 while 
the combined model’s average NMSE (without drift) is 1.0943 x 10~ 2 , 284 percent 
higher. Visually, however, y(n) and y,(n) for the combined model still match well. 
Therefore, the value of the dynamic feature of the combined model outweighs the 
disadvantage of the increased simulation error, and the combined model may still 
be used with reasonable confidence despite this increased error. 
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Figure 3.12: Fitted curves for inner pole pair, ‘42 with antenna simulator 
(10 pairs), (a) Magnitude, and (b) phase. 
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Figure 3.13: Second AR coefficient a 2 as a function of t and E n . 

Finally, the combined model can be related to the linear difference equa¬ 
tion with non-constant coefficients in Eq. 3.6. The simulation program stores only 
the roots of a (t,E n ) and b(£, E n ), not the coefficients themselves. However, these 
coefficients may easily be computed. For example, Fig. 3.13 shows <i 2 (t,E n ) for 
0 < t < 20 hours, and 2 x 10~ 6 < E n < 40 x 10 -6 watt-seconds. 

4. Adding the Dummy Load to the Combined Model 

Figures 3.14a and 3.14b show a data sequence pair for the '42 trans¬ 
mitter connected to the resistive dummy load instead of the antenna simulator. 
Figures 3.14c and 3.14d present sequences h(n) and hxs(n) and the pole-zero plot 
corresponding to hxn(n) for model order P = 2, Q = 1. The time-domain match is 
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Sample number, n 

(a) 


RF PULSE: OUMMY LQAO (PAIR 50) 



(b) 


Figure 3.14: ‘42 input and output, dummy load, pair 50. 
(b) RF pulse. 


(a) TDW and 
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POLE/ZERO PLOT (PAIR 50) 
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ORIGINAL & MOOEL SEQUENCES h(n) & hA8(n) 



Sompl* number, n 


(d) 


Figure 3.14: (c) Pole-zero plot of ‘42 LTI model, pair 50, with dummy 
load, and (d) original and model sequences, h(n) and ^(n). 



















not quite as good as with the antenna simulator sequences because a lower model 
order was chosen for the dummy load. Slightly higher model orders yield a modest 
decrease in MSE but the locations of the poles and zeros are no longer consistent 
and are thus unusable in the catenated model. As a result, order P = 2, Q = 1 is 
used and is considered adequate. Appendix C also includes curves for these dummy 
load poles and zeros. 

F. NONLINEAR, TIME-VARYING MODEL OF THE 
AN/FPN-44A TRANSMITTER 

The same procedure produced a combined model for the AN/FPN-44A trans¬ 
mitter, with order P = 6, Q = 5. A linear model of typical pair 71 has the following 
a and b: 

1.0000 
-5.9652 
14.8348 
a = -19.6972 

14.7246 
-5.8759 
0.9779 

The poles and zeros are 

‘ O.9947e j±0 - 0850 1 .OOOOeJ* 01023 ' 

poles = OMSW* 0 ** 59 zeros = 1.0025e J±00413 . (3.37) 

0.9956e> ±(M519 1.0255e'° 

All the data vectors for the ’44A were obtained on the same day, so typical time 
variations of this transmitter model were inferred. The combined model was tested 
with fifteen data pairs against the linear model from pair 71, as before. The linear 
model’s average NMSE is 1.3059 x 10 -J , and the combined model’s average NMSE 
is 5.0988 x 10 -3 . The catenated model thus reduces NMSE by 61.0 percent. 

Figures 3.15-3.17 are the plots for the ’44A modeling process. Appendix C 
contains the pole-zero curves for the ’44A. 
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Figure 3.15: Power spectrum of ‘44A with antenna simulator, pair 71. 
(a) TDW and (b) RF pulse. 
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Figure 3.15c: Closeup of power spectrum, ‘44A, pair 71. 
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Figure 3.16: (a) Pole-zero plot of ‘44A LTI model, pair 71, with antenna 
simulator, and (b) original and model sequences, h(n) and /ua(n). 
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POLE/ZERO PLOT (PAIR 85) 
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Figure 3.17: (c) Pole-zero plot of ‘44A LTI model, pair 85, with dummy 
load, and (d) original and model sequences h(n ) and 
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IV. THE STEEPEST DESCENT ALGORITHM 

A. INTRODUCTION 

This chapter presents am algorithm that generates and controls the Loran-C 
pulse shape automatically. First, the Loran-C control problem is described. Second, 
the algorithm is derived and its advantages and limitations are discussed. Finally, 
the algorithm is successfully tested with the combined ’42 and ’44A transmitter 
models from Chapter III. 

B. PULSE SHAPE CONTROL IN LORAN-C 

The VXIbus pulse shape control system (See Fig. 2.10) uses a batch pro¬ 
cessing approach to monitor and control the Loran-C pulse shape. The oscilloscope 
captures and stores in memory an entire RF pulse, or possibly an entire GRI or 
even a PCI. A resident computer program reduces the waveform(s) to parametric 
form, compares these parameters to the parameters from an ideal pulse, and pro¬ 
duces a new TDW in parametric form. In the final step, the program expands these 
TDW parameters into a digital data vector, which the AFG then converts into an 
analog TDW. This operation constitutes one iteration of the Loran-C pulse shape 
controller. The controller shapes the pulse by changing the TDW parameters to 
minimize the error in the RF pulse parameters. 

The simplest form of the controller uses a parametric form consisting of 16 
TDW half-cycle peak amplitudes 

*, = (4.1) 

and 16 RF pulse half-cycle peak amplitudes: 

(4.2) 
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These are tracked in successive control iterations; Fig. 4.1 illustrates this for x p (3) 
and y p (3). The half-cycle peak amplitudes of the ideal pulse of Eq. (2.1) are 

yop = (yop(l), yop(2), • • •, yo»(16)]', (4.3) 

with output error vector 

e y = y P -yop = [e y ( l),e v (2),---,e v (16)]'. (4.4) 

Thus the Loran-C pulse shape control problem is formulated as a 16-dimensional 
regulator problem where y p is maintained as close as possible to yop [Ref. 22]. 

C. THE STEEPEST DESCENT ALGORITHM 

1. Derivation 

A linear controller feedback which uses the method of steepest descent ef¬ 
fectively shapes the Loran-C pulse by minimizing the quadratic error e^We y , where 
W is a weighting matrix. In the steepest descent method, also used in adaptive fil¬ 
tering, the error gradient is used to find the bottom o f an error performance surface 
[Ref. 19, p. 197]. The controller presented here is appropriately called the steep¬ 
est descent algorithm in this thesis. This control algorithm, developed by Peterson 
and successfully implemented in a hardware simulation by Steinvorth, is derived as 
follows [Ref. 20, 21]. 

The transmitter is modeled as a linear system 

A 0 x„ = y p , (4.5) 

where A 0 is the matrix of impulse response samples (peak amplitude samples are 
used in this thesis): 

‘ fc p (0) 0 0 O' 

MD M°) o ••• o 

A 0 = ; . . . (4.6) 

. M16) MIS) MU) M°). 
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(b) 

Figure 4.1: Loran-C pulse shape control strategy, (a) Input and (b) out- 
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Let Xop be the vector of optimum TDW half-cycle peak amplitudes. Multiplying the 
true (but unknown) system matrix Ao and optimum TDW half-cycle peak vector 
Xop yields the ideal pulse peaks: 

A 0 Xop = yop • (4.7) 


Let 


®x — Xp Xop 


(4.8) 


be the vector of errors in the TDW half-cycle peak amplitudes. Beginning with A , 
the best estimate of A 0 , the first vector of TDW half-cycle peak amplitudes is 



(4.9) 


Vector Ax p updates x p at each iteration t in the direction of steepest descent on the 
quadratic error surface defined by e^We,: 


Ax, = -|v.,[e^We„), (4.10) 

where p is a small constant greater than zero and Ve z [«] is the vector of partial 
derivatives of scalar k with respect to the elements of e x . 

Equation (4.10) may be expanded to the form 


Ax, = -|v.,(elAfWA„eJ. (4.11) 

Diagonal matrix W weights the elements of e r and e„. W = I weights all the errors 
equally, while a W whose diagonal elements decrease in size from upper left to lower 
right implements a tighter tolerance in the first part of the pulse. Equation (4.11) 
then becomes 

Ax p = -pAg WA 0 e r , (4.12) 


then 


Ax p = -pAo We y . 


(4.13) 
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Matrix A, which is known, is substituted for unknown Ao, resulting in the steepest 
descent algorithm 

(4.14) 

For negatively phase coded pulses the ” sign is removed. Finally, if [A r WA 0 ] _1 
is positive definite and 


Ax p = -/iA T We y 


P < 


_ 2 _ 

largest eigenvalue of [A T WAo] ’ 


(4.15) 


the system is stable. 


2. Advantages and Limitations 


The algorithm of Eq. (4.14) has two main advantages and two pri¬ 
mary limitations. The algorithm’s advantages are its effectiveness and its simplicity. 
As the next section shows, the algorithm works well for the data available. Also, 
the simple parametric form allows the algorithm to be implemented with only two 
(16 x 16) matrix multiplications per iteration (with 768 real multiplications and 512 
real additions). The algorithm’s limitations are its assumption of an overall linear 
controller system and its inability to control zero-crossing times. By using small 
values of n and a standard 16-half-cycle TDW, the controller keeps the transmit¬ 
ter confined to an approximately linear region. If the controller is also “tuned” by 
updating A from time to time to compensate for time variations, the first limita¬ 
tion is not a serious one. The second limitation simply requires the zero-crossing 
times to be adjusted as they have always been: by keeping the transmitters well- 
maintained and well-tuned to 100 kHz. This limitation imposes no extra burden on 
Loran-C operation or maintenance, so likewise it is not a serious one. Admittedly, 
having an algorithm that adjusts zero-crossing times automatically would be quite 
advantageous, but this capability is not absolutely necessary in order to implement 
automatic pulse shaping successfully. 
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D. CONTROLLING PULSE SHAPE IN THE COMBINED 
MODEL USING THE STEEPEST DESCENT ALGO¬ 
RITHM 

The steepest descent algorithm effectively generates and controls the Loran-C 
pulse peak amplitudes in the combined models of both the ’42 and the ’44A. The 
zero-crossing tolerances were met for the ’42 but not the ’44A. 

Figures 4.2a through 4.2d illustrate the performance of the steepest descent 
algorithm with the ’42 combined model without time variations. Figure 4.2a shows 
the smooth convergence of the mean squared error e^e y for the first 100 iterations, 
and Fig. 4.2b presents the convergence of three measures of Loran-C error for the 
same 100 iterations. These are the half-cycle peak amplitudes (ensemble tolerance) 
and the half-cycle peak amplitudes (individual tolerance) for half-cycles 1-8 and 
for half-cycles 9-13. The temporary rise in mean squared error from iteration 10 
through 18 in Fig. 4.2a is a result of linearly controlling a nonlinear system. Figures 
4.2c and 4.2d show the first 800 samples of the ideal and synthetic RF pulses after 
100 iterations. All half-cycle peak amplitudes of the initial TDW (at iteration t = 1) 
were set to 0.4 volts. In these plots, normalized p 

largest eigenvalue of [A r WA] ,. , 

f*n = p -^ 1 - ( 4 - 16 ) 

is set to 0.7. Next, beginning with iteration 101 and ending with iteration 200, 
time variations are introduced at each iteration. This approximates a 25-hour pe¬ 
riod. These figures demonstrate the algorithm’s ability to compensate for slow time 
variations in the transmitter’s transfer function. 

Figures 4.4 and 4.5 illustrate similar success with the ’44A except that zero- 
crossing times are not within tolerance. Half-cycle peak amplitudes 15 and 16 of 
the final TDW are larger, to fill out the end of the controlled part of the pulse, but 
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Figure 4.2: Testing steepest descent algorithm with ‘42 and antenna, (a) 
Convergence of mean-squared error eje y , and (b) convergence of three 
measures of Loran-C error. 
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Figure 4.2: Testing steepest descent algorithm with ‘42 and antenna, (c) 
Final TDW and (d) ideal and synthetic RF pulses. 
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Iterations, t 


Figure 4.3c: Testing steepest descent algorithm with time-varying model, 
M2, drift parameters (cq coefficients). 

this has an unintended effect with the ’44 A combined model. This is because the 
tail drive decays exponentially starting from the amplitude of half-cycle 16. With 
such a large 16th half-cycle, the tail of the TDW contains almost as much energy as 
half-cycles 1-14. If actually implemented, half-cycles 15 and 16 could be attenuated 
to compensate for this problem. Also, the phase difference in Fig. 4.4d reflects 
the change in transmitter delays that may be produced by time variations in the 
transmitter’s transfer function. As stated before, emission delays are not addressed 
in this simulation. 

These figures show that the steepest descent algorithm effectively shapes the 
Loran-C pulse. However, these tests of the algorithm are unrealistic because com¬ 
plicating factors such as noise, quantization error and power supply droop are not 
included. Also, only one pulse of the PCI has been controlled. In the next chapter, 
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Magnitude 


the combined models of Chapter III and the steepest descent algorithm are incorpo¬ 
rated into a comprehensive simulation program which does provide a realistic test 
Experimental results from the simulation program are also included. 



Iteration*, t 

Figure 4.4a: Testing steepest descent algorithm with ‘44A and antenna. 
Convergence of mean-squared error e y e y . 
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Figure 4.4: Testing steepest-descent algorithm with 
(b) Convergence of three measures of Loran-C error 


‘44A and antenna, 
and (c) final TDW. 
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Figure 4.4d: Testing steepest descent algorithm with ‘44A and antenna. 
Ideal and synthetic RF pulses. 



Iteration*, t 

Figure 4.5a: Testing steepest descent algorithm with time-varying model, 
‘44A. Convergence of mean-squared error eje v . 
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(c) 

Figure 4.5: Testing steepest descent algorithm with time-varying model, 
‘44A. (b) Convergence of three measures of Loran-C error and (c) drift 
parameters (cq coefficients). 
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V. SIMULATION PROGRAM AND RESULTS 

A. INTRODUCTION 

In this chapter the models of Chapter III and the control algorithm of Chap¬ 
ter IV are incorporated into a comprehensive MATLAB computer program which 
simulates the pulse-shaping control process on the AN/FPN-42 and AN/FPN-44A 
transmitters. Also, key results obtained from this simulation program are featured. 

B. THE SIMULATION PROGRAM 

1. Structure 

The diagram in Fig. 5.1 shows the basic structure of the simulation 
program. In brief, the user selects the options for the simulation run, including 
which pulses of the PCI he or she wishes to control. For each of the selected 
pulses, the program completes a specified number of control iterations. A control 
iteration consists of obtaining the RF pulse, determining the error in the RF pulse 
parameters, and producing a new TDW. After the iterations are finished, a pulse 
analysis is performed and the program moves to the next selected pulse. This 
program simulates controlling the shorter rate of a dual-rated station, and from 
time to time the rate is blanked. When this occurs, the controller skips an entire 
control iteration and does not increment the loop counter. Thus the blanking process 
is simulated but is invisible to the pulse shape controller. 

2. Explanation of Features Appearing on Main Menu 
a. Main Menu 

The user controls the simulation program through a main menu, 
which appears in Fig. 5.2. Using this menu, the user configures the transmitter and 
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LORAN-C TUBE TRANSMITTER SIMULATION 
(c) Dean C. Bruckner & Murali Tummala, 1992 


2. Xmtr load 
4. Local ECD 
6. Imbalance 
8. Reset Xmtr 


Dummy Load 
0.00 us 
0.00 


1. Transmitter: AN/FPN-42 

3. Sampling freq: 10.00 MHz 
5. Resolution: 8 bits 

7. W. Noise pet: 0.39 

9. Pulses to control: [ 1 ] 

10. Pulses to analyze: [ 1 ] 

11. Number of iterations: 100 (1st pulse), 20 (following pulses) 

12. Xmtr parameter drift occurs every 0 iterations with norm mag. 1 

13. Xmtr switch occurs every 0 iterations (1st pulse only, when drift on) 

14. Display method: plot 

15. Control algorithm: Steepest Descent 

16. Display/change current algorithm parameters 

17. Access keyboard 18. Exit 

Enter number(s) to change parameters or <Enter> to begin: 


(e.g., l or [178]) 


Figure 5.2: Main menu, simulation program. 

control algorithm and selects the desired display, analysis and recording options. In 
this section each menu item is briefly explained. ' 

b. Transmitter Selection 

The user may select either the AN/FPN-42 or the AN/FPN-44A 
transmitter. The program loads the polynomial coefficients for the selected trans¬ 
mitter, which have been stored in a single matrix with one polynomial in each row, 
as in Fig. 5.3. The polynomial coefficients of the fcth root appear in adjacent rows — 
the first for magnitude and the second for phase. The program reinitializes variables 
governing the transmitter’s operation and resets the drive waveforms and control 
algorithm. 
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c. Transmitter Load 


The simulation program operates with either the antenna or the 
dummy load. The polynomial curves for the dummy load are in the lower partition 
of the matrix in Fig. 5.3. The program implements a load switch by resetting a row 
pointer for this matrix to select either the upper or lower partition. In its default 
mode, the program uses the dummy load to produce a near-optimum TDW for the 
antenna. The program switches to the antenna when the output errors fall below a 
threshold. This minimizes the time the pulse is out of tolerance when transmitting 
on the antenna. The “ideal” dummy load RF pulse used in this process was obtained 
by allowing the algorithm to converge on the antenna, switching to the dummy load 
and recording the output of this TDW. After switching to the antenna, usually the 
RF pulse is in tolerance within an iteration or two. Here the antenna and antenna 
simulator are used interchangeably. 


Row I | 

Pointer —► I I 

Magnitude c * c *-i 

Phase 1 

(root k) I I 




Antenna 


Dummy 

Load 


Figure 5.3: Polynomial coefficient matrix. 


d. Sampling Frequency 

This program runs at four data sampling frequencies: 10 MHz, 5 
MHz, 2.5 MHz, and 1.25 MHz, as discussed in Chapter III, Subsection D.l. The 
best error convergence is at f t = 10 MHz, but the program runs the fastest and 
requires the least storage at f, = 1.25 MHz. The algorithm resets the transmitter 
and control algorithm when a new sampling frequency is selected. 
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e. Local ECD 

The program controls the local (transmitted) ECD of the RF pulse 
by generating a new ideal Loran pulse with the desired ECD from Eq. (2.1) and 
using the new pulse in the control algorithm. Currently in the Coast Guard, ECD is 
controlled by inserting a phase shift called the Early Timing Adjust (ETA) into the 
TDW. This program bypasses the ETA altogether and successfully controls ECD 
to within 0.44/is in the range —2.5 ns < r < 2.5/iS by changing the ideal waveform. 
The LOIS program, used in the daily Loran-C system sample, is used to measure 
ECD by hand. 

f. Amplitude Resolution and System Noise 

The simulation program incorporates the noise model shown in Fig. 
5.4. The noise present in the actual data pairs may be duplicated in simulation by se¬ 
lecting eight-bit quantization and adding white noise to the synthetic TDW and 



Figure 5.4: Transmitter system noise model. 
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TABLE 5.1: AVERAGE SNR OF MEASURED DATA PAIRS 
(ANTENNA) 


Transmitter 

SNR (dB) 

TDW 

RF Pulse 

AN/FPN-42 

56.4 

62.0 

AN/FPN-44A 

56.1 

66.6 


RF pulse until the SNRs of both match the average SNRs of the actual data (as 
defined in Chapter III, Subsection B.2.b.). These average SNRs are listed in Table 
5.1. 

Because the relative amplitudes of the ’42 and ’44A waveforms are 
different, the standard deviation of the white noise is expressed as a percentage 
of the maximum positive amplitude of the waveform. The SNRs in Table 5.1 are 
achieved in simulation using the settings in Table 5.2. 

The user specifies the number of bits and the noise percentage of 
the RF pulse in menu items five and seven, respectively; the program then sets the 
TDW noise percentage automatically by multiplying the RF pulse noise percentage 
by 2.7 for the ‘42 and 1.8 for the ‘44A. Other quantization settings available are 
S q = 12 bits, 5, = 16 bits and S q — oo (maximum resolution, to machine precision). 
These represent a best-case scenario, because all the quantization levels are used. 
In the real system, some of the levels at the top and bottom are not usually used to 
avoid saturation, reducing the effective bit resolution. The capability to reproduce 
the observed noise level in the transmitter system is extremely important as it allows 
the simulation to be a realistic one. Results of different quantization settings are 
presented in Section C of this chapter. 






TABLE 5.2: PROGRAM SETTINGS WHICH REPRODUCE SNR OF 
MEASURED DATA 


Transmitter 

Number of 
bits, S q 

Std. Dev. of White Noise 
(% of peak amplitude) 

TDW 

RF Pulse 

AN/FPN-42 

8 

1.05 

0.39 

AN/FPN-44A 

8 

0.97 

0.54 


g. Transmitter Imbalance 

As described in Subsection C.l.e. of Chapter II, an imbalance be¬ 
tween the two vacuum-tube amplifier banks in a transmitter reduces the amplitude 
of the negatively phase coded pulses. The program simulates this imbalance by 
reducing the amplitudes of these RF pulses by a percentage defined by the user 
in this menu item. The program automatically compensates for this imbalance by 
increasing the TDW amplitude by an appropriate amount. As with ETA, the phase 
code balance adjustment in the PGEN is bypassed entirely. 

h. Reset Transmitter 

When the program completes controlling and analyzing all the se¬ 
lected pulses in a PCI, the main menu appears again and the user has the option to 
continue where the program left off. The reset feature allows the user to start the 
control process from the beginning again without exiting the program. When the 
user selects this item, the program resets the drive waveforms, the control algorithm 
and the random transmitter drift settings (if drift is enabled), but it leaves intact 
the other control and analysis settings in the main menu. 
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i. Pulses to Control 

The program can control any or all of the pulses in a PCI, as speci¬ 
fied by the user. The TDWs for all selected pulses of the PCI are stored in successive 
columns of matrix D c , which represents a data output buffer to the AFG. The con¬ 
trol approach is sequential, beginning with the first selected pulse. The program 
“drives” the transmitter model by presenting the TDW in column one of D c as an 
input argument to the function XMTR. The resulting RF pulse is in turn presented 
as an input argument to the control algorithm, which produces the new TDW. This 
TDW, which is the best estimate of the optimal TDW for each pulse, is loaded into 
all the columns of D c and proper phase-coding is applied. The amplitude of each 
TDW may also be scaled up exponentially to compensate for power supply droop as 
explained later in this section. When the specified number of control iterations are 
completed, the final RF pulse is stored in column one of matrix R c and the program 
moves to the next selected pulse. As the program controls the pth pulse, columns p 
and following of D c are updated every iteration, but columns one through p — 1 are 
not. When the entire process is completed, matrix D c contains the best estimates 
of the optimal TDWs for all selected pulses, and R c contains the selected output 
RF pulses. In the VXIbus system, the output data buffer can easily be dumped to 
the AFG. With an MPT to set the proper time of emission, the desired TDW would 
be sent to the transmitter. 

j. Pulses to Analyze 

When the program finishes controlling an RF pulse, it performs an 
analysis of that pulse and the control process that produced it. The program’s 
default setting is to analyze every selected pulse, but the user may suppress any 
or all of these analyses. The program then prints the results to a screen and to 
an ASCII text file, as in Fig. 5.5. Next, the program plots the Loran-C errors and 
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LORAN-C PULSE ANALYSIS 


System Parameters: 


1. Transmitter: 

AN/FPN-42 

2. Xmtr load: 

Antenna 

3. Sampling freq: 

10.00 MHz 

4. Local ECD: 

0.00 us 

5. Resolution: 

0 bits 

6. Imbalance: 

0.00 

7. W. Noise pet: 

0.00 




8. Pulse 1 

9. Number of iterations: 100 

10. Xmtr parameter drift occurred every 0 iterations w/ norm mag. l 

11. Xmtr switch occurred every 0 iterations 

Parameters for control algorithm: Steepest Descent 

1. Initial Mu (duosny load) : 0.8 

2. Mu after load switch (antenna): 0.7 

3. Mu_max: 0.0008492 

4. Weighting Matrix: W = I 
Press <Enter> to continue 


PULSE IN TOLERANCE (BCD & power spectrum not checked) 

Press <Enter> to continue 

MSB out / Ens err / MaxE 1-8 / MaxE 9-13 
ans » 

0.0053 0.0035 0.0083 0.0050 

errjmean » 

0.0057 0.0035 0.0083 0.0047 

err_sdev « 

2.1403e-04 1.7029e-06 3.2212e-05 1.9595e-06 

Peak amplitudes in tolerance for all iterations examined 

Ratel blanked before the following iteration numbers: 
blanksav . 6 20 27 33 40 54 63 77 91 

Avg time per iteration: 1.655 seconds 

Press <Bnter> to continue 


Figure 5.5: Pulse analysis printout. 
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Normalized pulse peak values 

average, last iter, ideal, diff (1-3), diff (2-3) 


ans 3 


0.0074 

0.0074 

0.0157 

-0.0083 

-0.0083 

-0.0793 

-0.0793 

-0.0833 

0.0041 

0.0041 

0.1915 

0.1915 

0.1901 

0.0013 

0.0014 

-0.3180 

-0.3181 

-0.3158 

- 0.0022 

-0.0023 

0.4470 

0.4471 

0.4454 

0.0016 

0.0017 

-0.5711 

-0.5712 

-0.5696 

-0.0015 

-0.0016 

0.6828 

0.6828 

0.6813 

0.0014 

0.0015 

-0.7782 

-0.7781 

-0.7771 

- 0.0010 

- 0.0010 

0.8584 

0.8584 

0.8556 

0.0028 

0.0028 

-0.9211 

-0.9214 

-0.9164 

-0.0047 

-0.0050 

0.9625 

0.9631 

0.9598 

0.0027 

0.0033 

-0.9856 

-0.9861 

-0.9872 

0.0015 

0.0011 

1.0000 

1.0000 

1.0000 

0 

0 

-1.0067 

-1.0064 

- 1.0001 

-0.0066 

-0.0063 

0.9965 

0.9965 

0.9892 

0.0073 

0.0073 

-0.9668 

-0.9675 

-0.9692 

0.0024 

0.0017 


SNR (tdw) = 83.93 dB 
SNR (rf) = 108.3 dB 


Figure 5.5 (continued): Pulse analysis printout. 

output mean squared error, the final TDW and RF pulse, the rest times and esti¬ 
mated power supply voltage (explained later in this section), and the drift parame¬ 
ters. These plots are also recorded in META file format, 
k. Number of Iterations 

Here the user selects the number of control iterations for the first 
selected pulse and for all following pulses. The default setting for the first selected 
pulse is 100, which is usually sufficient time for the algorithm to converge. Because 
the TDWs of the following pulses are also replaced every iteration during the con¬ 
vergence of the first selected pulse, the pulses are in tolerance or nearly in tolerance 
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at their first control iteration. Bringing these errors into convergence thus requires 
fewer iterations and so the default setting is twenty iterations. 

l. Transmitter Drift 

In this item, the user sets both the time scale and the magnitude 
of the transmitter parameter drift described in Chapter III, Subsections E.2. and 
E.3. Setting tj = 300 and the normalized magnitude equal to one reasonably ap¬ 
proximates the time variations in an actual transmitter. Reducing the value of or 
increasing the normalized drift magnitude artificially increases these time variations, 
which is useful in testing and analysis. 

m. Transmitter Switch 

By introducing a random step change into the cq parameters shown 
in Fig. 5.3, a switch to a different transmitter may be simulated. The program 
performs transmitter switches at a regular interval, as often as the user desires. In 
this way, the algorithm may be tested with an infinite number of different ’42 and 
’44A transmitter models as no two random settings of the cq parameters are exactly 
alike. When the number of iterations per transmitter switch is set to zero, however, 
the identical transmitter model is produced every time. To obtain a random setting 
and to disable time variations, the user sets the number of iterations per transmitter 
switch to the number of control iterations (menu item 11) plus one. The load remains 
the same during a transmitter switch. The performance of the control algorithm 
following transmitter switches appears in Section C of this chapter. 

n. Display Method 

To monitor the errors as the algorithm converges, the user may 
select the text line option of Fig. 5.6a or the plot option of Fig. 5.6b. Or, the 
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Iter # 

/ MSB out 

/ Bns err 

/ MaxE 1-8 

/ MaxE 9-13 

top 

a 







I.0000 

0.4163 


0.0219 

0.0356 

0.0230 


Iter # 

/ MSB out 

/ 

Bns err 

/ MaxE 1-8 

/ MaxE 9-13 

top 

a 







2.0000 

0.4163 


0.0219 

0.0356 

0.0230 


Iter # 

/ MSB out 

/ 

Bns err 

/ MaxE 1-8 

/ MaxE 9-13 

top 

a 







3.0000 

0.3029 


0.0268 

0.0428 

0.0160 


Iter # 

/ MSB out 

/ 

Bns err 

/ MaxE 1-8 

/ MaxE 9-13 

top 

m 







4.0000 

0.2578 


0.0266 

0.0434 

0.0128 


Conv#rg«nc* of En*«mbl« Error, Me* «rra 1-8, Mo* errs 9-13 




nioi»««05i ! bits*0 


50 60 70 80 9 

lt«rotion», t 

(b) 

Figure 5.6: Options for displaying error convergence, (a) Text line, and 
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STEEPEST DESCENT ALGORITHM PARAMETERS 
(Written by CAPT B. B. Peterson, 
Modified by Deem Bruckner & Murali Tummala) 


1. Initial Mu (dummy load): 0.8 

2. Mu after load switch (antenna): 0.3 

3. Weighting Matrix: W = I 

4. Access keyboard 

Enter number(s) to change parameters or <Enter> to return.- 
(e.g., 1 or [134]) ====> 


Figure 5.7: Menu for steepest descent algorithm. 

user may suppress both of these. The plot is the most comprehensive display 
method. Load switches are marked by an “a” and transmitter switches by an “z”. 

o. Control Algorithm 

Currently only the steepest descent algorithm has been implemented 
in this simulation program. Because the program is constructed in modular form, 
other algorithms may easily be added. In this case, the user could select a new 
algorithm from the main menu. 

p. Display/Change Current Parameters 

The parameters for each algorithm are changed using a submenu. 
Figure 5.7 displays this menu for the steepest descent algorithm. 
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q. Access Keyboard 

In this option, the user accesses the MATLAB command line. Here 
he or she may examine or change parameters not accessible at the menus. 

3. Explanation of Features Not Appearing in the Main Menu 
a. Simulating Power Supply Droop 

The problem of power supply droop has a negative effect on a closed 
loop Loran-C control system, as described in Subsection B.2.d. of Chapter II. The 
amplitude of each successive pulse in a GRI is reduced by a small amount, because 
the power supply of the transmitter does not fully recover within the 1000 /xs pulse- 
to-pulse interval. This program models the power supply voltage v v , as a charging 
exponential with four parameters: 

• The normalized power supply voltage just prior to the last pulse (v t , in the 
range (0,1]), 

• The decrease in power supply voltage due to transmitting the last pulse (Au/), 

• The time the power supply has rested since the last pulse ( t v ), and 

• The time constant of the charging exponential (r„). 

These are shown graphically in Fig. 5.8. It was first assumed that the transfer 
function of the transmitter changed measurably with decreases in the power supply 
voltage; however, an analysis of an entire GRI revealed no predictable changes in 
pole and zero locations. Therefore, the effect of droop is modeled solely as a decrease 
in RF pulse amplitude. 

The quantity v pt is estimated before obtaining each RF pulse. In 
the function XMTR, the RF pulse is multiplied directly by the v pe estimated for 
that pulse. This simulates the decrease in amplitude due to droop. Quantities At>* 
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R«eov«ry of pp-*<olt for AN/PPN-42 



Figure 5.8: Charging Exponential model of power supply voltage. 


and r v are constant for the entire simulation, but quantities v t and t v are not. Here 
a simplification is made: the power supply voltages of pulses two through eight in 
each GRI are assumed to decrease exponentially from pulse one, according to the 
relation 




= 6>V n 


puke p 


puke 1 


P = 2,---,8, 


(5.1) 


where 6 = 0.992 for the ’42 and 6 = 0.9988 for the ’44A. Now v pi must be estimated 


for the first pulse of every control iteration in order to find v p , for any of the other 
pulses. 


Because pulse eight of the last GRI transmitted always precedes 
pulse one, vt is held constant in the simulation at 0.935 for the '42 and 0.99 for the 
? 44A. This reflects measured values of droop of 6.5 percent for the ’42 and 1.0 percent 


109 








for the ’44A. Realistic values of variable t v for each control iteration, however, are 
found only by simulating dual-rated operation. 

b. Simulating Dual-Rated Operation 

As explained in Chapter II, Subsection B.4., a dual-rated station 
transmits GRIs for two adjacent chains, and the two GRIs are not synchronized 
to each other. When the blanking intervals enclosing the GRIs overlap, one of the 
rates is blanked. When the blanking intervals are very close to each other but do 
not overlap, however, the droop of the power supply is accentuated. In this case, 
the transmitter transmits eight pulses spaced 1000 ps apart followed by a rest as 
short as 1900 ps, which is then followed by eight more pulses spaced at 1000 /is. 
The interval betw r een the two GRIs is the quantity t v . 

A problem arises when the VXIbus control system happens to sam¬ 
ple a pulse in this second GRI. The peak amplitude of the sampled pulse from this 
GRI will be noticeably smaller than the average peak amplitude of that pulse from 
other GRIs. This introduces a transient disturbance into the control algorithm. The 
purpose of including dual-rate operation in this simulation is to study the response 
of the control algorithm to these transient effects. 

The MATLAB function REST was written to produce a vector of 
consecutive rest times t v , in seconds, prior to the GRIs of the shorter rate in a 
dual-rated station. For two secondary rates 99,400 ps and 59,900 ps, 200 samples 
of this vector are plotted in Fig. 5.9. From this figure, the periodic nature of this 
function is clear. In this program the pulse shape controller captures a pulse from 
a GRI every four seconds; with 66 GRIs in a four-second period, realistic values of 
f v for pulse one in each GRI are thus obtained by taking every 66th sample from 
the vector of rest times. Another vector of the same length b v indicates the samples 
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CONSECUTIVE REST TIMES 



Group Repetition Intervals 

Figure 5.9: Vector of consecutive rest times, 200 samples. 

for which the shorter rate was blanked; if the controller encounters one of these, it 
records the blanked iteration (see Fig. 5.5) and skips that GRI altogether. 

Vector t„ is 503 samples long; when the end of the vector is reached 
the counter wraps around to the beginning again. Because 503 is a prime number, 
all of the samples in the vector should eventually be used if the simulation continues 
long enough. Figures 5.10a and 5.10b show the error convergence plots for a ’42 
simulation run of 100 samples. The controller switches from the dummy load to 
the antenna at iteration 24. Figures 5.10c and 5.10d show the estimated values 
of t v and v pt for each iteration. The increase in error is clearly visible between 
iterations 29-32. Fortunately, the algorithm is able to compensate and these errors 
are transient. 
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Convergence of Ensemble Error, Max errs 1-8, Mcx errs 9-13 



0 10 20 30 40 50 60 70 80 90 100 
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Figure 5.10: Effect of transient disturbances, (a) Mean squared error and 
(b) three measures of Loran error. 
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Figure 5.10: Effect of transient disturbances, continued, (c) Rest times 
and (d) estimated power supply voltage. 






































c. Exponential Scaling of TDW Amplitudes 

In order to speed the convergence of pulses two through eight of the 
GRI, the program anticipates the droop of the transmitter. It boosts the ampli¬ 
tude of each successive TDW slightly when storing the TDWs in buffer matrix D c 
according to the relation 


(TDW) I = P P (TDW) 


P = 


(5.2) 


Ipulse p I puke 1 

where — 1.02 for the ’42 and = 1.0015 for the ’44A. These values are not exactly 
the inverses of the 6 values given earlier; no specific a priori knowledge of droop was 
assumed, and these values were obtained experimentally. This feature does in fact 
speed convergence; in most cases, pulses two through eight converge within two to 
three iterations. 


C. RESULTS 

X. Realistic Simulation of the AN/FPN-42 

In a realistic simulation of the ’42 with 600 control iterations, the steepest 
descent algorithr shaped the Loran pulse effectively but maintained the pulse peaks 
in tolerance for only 78.3 percent of the control iterations in a test interval beginning 
at interation 101 and ending at iteration 600. The zero-crossings and ECD of the 
final RF pulse were in tolerance, however. The system noise drove the pulse peaks 
out of tolerance repeatedly. These out of tolerance cases would not necessarily 
require blink, but they are undesirable. Whether or not this would be acceptable 
in actual Loran operation is a policy matter for the Coast Guard to decide. In 
any case, these fluctuations are an indication that the steepest descent algorithm is 
sensitive to noise and that it is not as robust as perhaps it should be. Figures 5.11a 
through 5.1 lg are the plots for the first 100 iterations of this run; Fig. 5.1 lh is the 
printout for the test interval. 
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Figure 5.11: Realistic simulation results, ‘42. (a) Mean squared error and 
fbT three measures of Loran error. 
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Figure 5.11: Realistic simulation results, ‘42. (c) Final TDW and (d) 
ideal and final synthetic RF pulses. 
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Figure 5.11: Realistic simulation results, ‘42. (e) Rest times and (f) 

estimated power supply voltages. 
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LORAN-C PULSE ANALYSIS 


System Parameters: 


1. Transmit ter: 

AN/FPN-42 

2. Xmtr load: 

Antenna 

3. Sampling freq: 

10.00 MHz 

4. Local ECD: 

0.00 us 

5. Resolution: 

8 bits 

6. Imbalance: 

0.00 

7. W. Noise pet: 

0.39 




8. Pulse 1 

9. Number of iterations: 100 

10. Xmtr parameter drift occurred every 0 iterations w/ norm mag. 1 

11. Xmtr switch occurred every 0 iterations 

Parameters for control algorithm: Steepest Descent 

1. Initial Mu (dummy load): 0.8 

2. Mu after load switch (antenna): 0.3 

3. Mu_max: 0.0008492 

4. Weighting Matrix: W = I 
Press <Bnter> to continue 


PULSE IN TOLERANCE (BCD & power spectrum not checked) 

Press <Enter> to continue 

MSB out / Ens err / MaxE 1-8 / MaxE 9-13 
ans « 

0.0212 0.0047 0.0085 0.0133 

err mean « 

“0.0249 0.0069 0.0113 0.0092 

err sdev * 

T.8736e-02 3.3487e-02 4.3991e-03 4.0186e-03 

Peak anqplitudes in tolerance for 83.6 % of iterations examined 

Ratel blanked before the following iteration numbers: 
blanksav > 

10 23 37 51 56 70 81 83 94 

Avg time per iteration: 3.59 seconds 

Press <Enter> to continue 


Figure 5.11h: Realistic simulation results,‘42, pulse analysis printout. 
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Normalized pulse peak values 


average. 

last iter. 

ideal, 

diff (1-3), 

diff (2 

0.0127 

0.0156 

0.0157 

-0.0030 

- 0.0000 

-0.0829 

-0.0859 

-0.0833 

0.0005 

-0.0026 

0.1932 

0.1953 

0.1901 

0.0031 

0.0052 

-0.3192 

-0.3203 

•0.3158 

-0.0034 

-0.0045 

0.4473 

0.4453 

0.4454 

0.0018 

-0.0001 

-0.5708 

-0.5781 

-0.5696 

-0.0012 

-0.0085 

0.6827 

0.6875 

0.6813 

0.0013 

0.0062 

-0.7778 

-0.7734 

-0.7771 

-0.0007 

0.0037 

0.8584 

0.8672 

0.8556 

0.0028 

0.0116 

-0.9225 

-0.9297 

-0.9164 

-0.0061 

-0.0133 

0.9637 

0.9688 

0.9598 

0.0039 

0.0089 

-0.9863 

-0.9844 

-0.9872 

0.0008 

0.0028 

1.0000 

1.0000 

1.0000 

0 

0 

-1.0084 

-1.0078 

-1.0001 

-0.0083 

-0.0077 

0.9979 

0.9922 

0.9892 

0.0087 

0.0030 

-0.9657 

-0.9609 

-0.9692 

0.0034 

0.0082 


SNR (tdw) * 54.54 dB 
SNR (rf) > 62.58 dB 


Figure 5.11h (continued): Pulse analysis printout. 

2. Realistic Simulation of the AN/FPN-44A 

In a similar realistic simulation of the ’44A, the steepest descent algo¬ 
rithm effectively shaped the Loran pulse and maintained the pulse peaks in tolerance 
for 96.5 percent of the control iterations in the test interval. The ECD of the final 
RF pulse is in tolerance but the second zero-crossing of this pulse is not in tolerance. 
Figures 5.12a through 5.12h contain the plots and printout of this simulation run, 
as before. 

3. Controlling ECD 

The steepest descent algorithm effectively controlled the ECD of the av¬ 
erage half-cycle peak amplitudes to within 0.44 /is, as Table 5.3 shows, during a 
test interval beginning 20 samples after the antenna switch and ending with sample 
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Figure 5.12: Realistic simulation results, ‘44A. (a) Mean squared error 
and (b) three measures of Loran error. 
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Figure 5.12: Realistic simulation results, ‘44A. (c) Final TDW and (d) 
ideal and final synthetic RF pulses. 
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Figure 5.12: Realistic simulation results, ‘44A. (e) Rest times and (f) 
estimated power supply voltages. 
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LORAN-C PULSE ANALYSIS 


System Parameters: 


1. Transmitter: 

AN/FPN-44A 

2. Xmtr load: 

Antenna 

3. Sampling freq: 

10.00 MHz 

4. Local BCD: 

0.00 us 

5. Resolution: 

8 bits 

6. Imbalance: 

0.00 

7. W. Noise pet: 

0.54 




8. Pulse 1 

9. Number of iterations: 100 

10. Xmtr parameter drift occurred every 0 iterations w/ norm mag. l 

11. Xmtr switch occurred every 0 iterations 

Parameters for control algorithm: Steepest Descent 

1. Initial Mu (dummy load): 0.8 

2. Mu after load switch (antenna): 0.3 

3. Mu_max: 0.217 

4. Weighting Matrix: W • I 
Press <Bnter> to continue 

Zero crossings exceed limits by following amounts (ns): 
z_err * 

-« 

-44.0257 
-37.4828 
0 
0 
0 
0 
0 
0 
0 
0 
0 

PULSE OUT OF TOLERANCE 
Press <Enter> to continue 

MSB out / Ens err / MaxE 1-8 / MaxE 9-13 

ans « 

0.0010 0.0048 0.0104 0.0041 

err mean = 

”0.0014 0.0074 0.0138 0.0108 

err sdev ■ 

3.0317e-04 1.7569e-03 3.7S97e-03 4.473Ge-03 

Peak amplitudes in tolerance for 93 % of iterations examined 

Ratel blanked before the following iteration numbers: 
blanksav • 

8 11 21 24 25 39 50 52 63 77 82 

Avg time per iteration: 2.496 seconds 

Press <Enter> to continue 


Figure 5.12h: Realistic simulation results, ‘44A, pulse analysis printout 
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Normalized pulse peak values 
average, last iter, ideal, 

diff (1-3), 

diff (2 

ans a 

0.0280 

0.0261 

0.0157 

0.0123 

0.0104 

-0.0781 

-0.0783 

-0.0833 

0.0053 

0.0051 

0.1861 

0.1913 

0.1901 

-0.0040 

0.0012 

-0.3178 

-0.3130 

-0.3158 

-0.0020 

0.0028 

0.4458 

0.4435 

0.4454 

0.0004 

-0.0020 

-0.5639 

-0.5652 

-0.5696 

0.0057 

0.0044 

0.6761 

0.6783 

0.6813 

-0.0052 

-0.0031 

-0.7749 

-0.7739 

-0.7771 

0.0022 

0.0032 

0.8534 

0.8522 

0.8556 

-0.0023 

-0.0034 

-0.9099 

-0.9130 

-0.9164 

0.0065 

0.0033 

0.9533 

0.9565 

0.9598 

-0.0065 

-0.0033 

-0.9855 

-0.9913 

-0.9872 

0.0017 

-0.0041 

1.0000 

1.0000 

1.0000 

0 

0 

-0.9899 

-0.9913 

-1.0001 

0.0101 

0.0088 

0.9753 

0.9739 

0.9892 

-0.0139 

-0.0153 

-0.9763 

-0.9826 

-0.9692 

-0.0071 

-0.0134 


SNR (tdw) = 56.76 dB 
SNR (rf) = 66.84 dB 


Figure 5.12h (continued): Pulse analysis printout. 

Table 5.3: CONTROLLING ECD (ALL VALUES IN fxs, WITH 8-BIT 
RESOLUTION AND WHITE NOISE ADDED) 





Measured 



Desired 

ECD of 

ECD 

Transmitter 

ECD 

Avg. y 

Error 


AN/FPN-42 



AN/FPN-44A -2.50 

0 

2.50 


-2.32 

0.04 

2.20 


-2.06 

0.01 

2.42 



.44 
.01 
-0.08 
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100. The realistic settings of Table 5.2 were used here also, with = 0.7. ECD 
values were measured by hand using the LOIS program, which accepts the first 
eight half-cycle amplitudes as its sole input arguments. The error converged to 
lower values for positive ECDs than for negative ECDs. 

Because ECD was not computed for each iteration, the maximum ECD 
error and the variance of the ECD errors are not available. Therefore, the precise 
effects of noise and quantization error on ECD are not known. 

4. Performance Improvement With Greater Bit Resolution 

Table 5.4 presents the improved performance with the ’42 that results 
from adding more quantization levels to the 5,-bit quantizer in the noise model of 
Fig. 5.4. A 600-iteration test was conducted as in Section B of this chapter, with 
the same test interval (iterations 101-600). Bit resolutions S q = 7 and S q = 10 
were set manually using the keyboard. Resulution S q = 7 represents a worst case 
for the VXIbus system, where the waveform's amplitude spans only half the vertical 
oscilloscope scale and uses only 128 of 256 quantization levels available. Resolution 
S q = 8 represents the best case where all 256 levels are used. The actual performance 
of the system should lie in between these two, closer to S q = 8. 

Three measures of the algorithm’s performance were selected for this 
comparison: the SNRs of both TDW and RF pulse, averaged over all iterations 
in the test interval; the percentage of iterations in which the pulse peaks are in 
tolerance; and the mean and standard deviation of the half-cycle amplitude error 
(ensemble tolerance). The maximum allowable value for this third measure is 0.01. 
The ensemble tolerance was chosen because it is usually exceeded first out of all of 
the pulse peak amplitude tolerances. The variance of the white noise remains the 
same, from Table 5.2. 
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Table 5.4: PERFORMANCE IMPROVEMENT WITH GREATER BIT 
RESOLUTION (‘42 TRANSMITTER, WITH WHITE NOISE ADDED 
AND fi n = 0.3) 


Bits 

S, 

Average 

SNR 

%of 
Iter, 
in Tol. 
(peaks) 

Half-Cycle Peak 
Ens. Error 

TDW 

RF 

Mean 

Std. Dev. 

7 

54.0 

59.9 

79.8 

0.0085 

0.0036 

8 

54.2 

62.1 

78.3 

0.0077 

0.0037 

10 

54.2 

63.3 

87.1 

0.0065 

0.0039 

12 

54.2 

63.4 

86.3 

0.0065 

0.0034 

16 

54.2 

63.4 

81.9 

0.0070 

0.0037 

oo* 

54.2 

63.4 

81.2 

0.0071 

0.0040 


* To machine precision. 

Although the accuracy of this simulation where S q / 8 is not supported 
by any data (to date, no oscilloscope for the VXIbus exists with more than eight 
bits), Table 5.4 shows a marginal improvement in all measures with S q = 10, but 
the performance worsened for higher values of S q . This is because quantizing into 
discrete levels removes completely the noise with amplitude less than one-half the 
bin width. Additional filtering is thus appropriate for all values of S q . Similar results 
should be expected for the ‘44A. 

5. Performance Improvement with Less Noise 

Reducing the variance of the white noise allowed the steepest descent 
algorithm to keep the pulse peaks in tolerance consistently with the ’42, while still 
using quantization with 8 bits. Using the same test procedure (used in Chapter IV 
and in the first part of this chapter) reducing fi n from 0.7 to 0.3 did reduce the 
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steady-state errors but did not increase the SNRs. The improved performance is 
apparent in Table 5.5. 


TABLE 5.5: PERFORMANCE IMPROVEMENT WITH LESS NOISE 
(‘42 TRANSMITTER, WITH 8-BIT RESOLUTION AND p n = 0.3) 


RF 

Noise 

(%) 

Average 

SNR 

% of Iter, 
in Tol. 
(peaks) 

Half-Cycle 

Peak Ens. Error 

TDW 

RF 

Mean 

Std. Dev. 

0.39 

54.1 

62.1 

78.3 

0 .0077 

0.0037 

0.25 

57.9 

64.7 

92.5 

0.0060 

0.0024 

0.20 

59.8 

65.7 

95.2 

0.0056 

0.0020 

0.15 

62.1 

66.7 

98.3 

0.0053 

0.0014 

0.10 

65.3 

67.9 

100 

0.0049 

0.0009 


* To precision of machine 


An acceptable solution may not necessarily require the pulse peaks to 
meet the criteria in the signal specification for 100 percent of control iterations. The 
minimum acceptable performance level and the method of reducing system noise are 
thus both subjects for further study. A similar performance improvement may be 
expected with the ’44A. 

6. Behavior Following a Transmitter Switch 

The steepest descent algorithm converges to the same degree as above 
after switching to a different ’42 transmitter of which it has no specific a priori 
knowledge. In Fig. 5.13, iterations 1-75 show the error convergence for a trans¬ 
mitter with realistic settings whose Co coefficients have been randomly displaced 
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as discussed in Subsection D.2. of Chapter III. At iteration 76, the coefficients are 
changed again, simulating a switch to a different transmitter. The algorithm’s quick 
response shows that it is versatile enough to use with transmitters at other Loran 
stations in the Coast Guard. Similar performance was seen with the ’44A. Again, 
system noise prevents the error from converging any further. 

7. Controlling the Entire PCI 

The steepest descent algorithm also performed well in controlling an 
entire PCI in a realistic simulation with a 3.0 percent transmitter imbalance intro¬ 
duced as well. Here /x n = 0.7. The pulse-to-pulse ECD and amplitude tolerances 
from Chapter II, Subsection B.2.e. were easily met. The shortcomings of the con¬ 
trol algorithm in dealing with noise appear also in pulses two and following, but the 
features of the program designed to control these pulses worked successfully. The 
program controlled ECD without ETA, compensated for droop without a separate 
droop correction circuit, and corrected phase code bounce without a separate phase 
code balance adjustment. 
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Mogratude 


Convorgenco of E«*«mbl« Error, Mox «rr* 1-8, Mox orr* 9-13 



Figure 5.13: Algorithm performance after transmitter switch at 

iteration 76. 


131 


















[THIS PAGE INTENTIONALLY LEFT BLANK] 


132 





VI. CONCLUSIONS 


A. CONCLUSIONS 

Modernizing the control systems for Loran-C vacuum-tube transmitters re¬ 
quires a control algorithm to shape the Loran pulse automatically, and in this thesis 
an algorithm was adapted for this purpose. In order to test the algorithm fully and 
to provide a tool for future study, a detailed simulation program for two classes of 
tube transmitters, the AN/FPN-42 and AN/FPN-44A, was developed. This pro¬ 
gram incorporates discrete-time HR models of each transmitter. 

Based on an initial assumption of LTI behavior at a given operating point, 
a linear difference equation with non-constant coefficients was chosen to represent 
the dynamics of the transmitters. Frequency-domain deconvolution, in conjunction 
with median smoothing, yielded an accurate estimate of the unit sample response 
at each operating point. Next, the least squares modified Yule-Walker method and 
Shank’s method provided a non-minimum phase pole-zero model of each sequence. 
These models were catenated to represent the transmitter’s nonlinearities, and time 
variations were added to form a combined model. The non-constant coefficients 
of the difference equation were defined as functions of time and the energy of the 
normalized TDW. The accuracy of this model was then demonstrated for both the 
’42 and the ’44A transmitters. 

Next, a linear feedback controller which uses the method of steepest descent to 
minimize the quadratic output error was derived, and its advantages and limitations 
were discussed. The algorithm successfully shaped the pulse with both the ’42 and 
the ’44A by bringing the pulse peaks into tolerance, although zero-crossing tolerances 
were exceeded in some cases. Then, the models and the control algorithm were 
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incorporated into a simulation program. To this program were added the details of 
the Loran transmitter system which affect pulse shape. Finally, the algorithm was 
tested in a variety of transmitter system settings and behaviors. From these tests, 
four main conclusions can be made. 

First, based on all the data available, the MATLAB simulation program and 
the nonlinear, time-varying models it contains accurately represent the behavior of 
the ’42 and ’44A transmitter systems over the range of operating points used. The 
assumption of LTI behavior at each operating point is a valid one, and the model 
reproduces it faithfully. The details of Loran operation added to the program make 
the simulation a realistic one. Therefore the results obtained are directly applicable 
to the VXIbus system. 

Second, the steepest descent algorithm shapes the pulse effectively in realistic 
simulations of both the ’42 and ’44A transmitters, with two significant shortcomings: 
the zero-crossing tolerances are exceeded occasionally with the ’42 and always with 
the ’44A, and the algorithm is sensitive to system noise. This noise drives the pulse 
peaks out of tolerance frequently. Still, under these conditions, the algorithm kept 
the ECD of pulse one in tolerance and quickly produced an entire PCI which met the 
pulse group tolerances for amplitude and ECD, even when a transmitter imbalance 
was added. Further, the algorithm reshaped the pulse effectively after transmitter 
switches. 

Neither of these two shortcomings necessarily disqualify the algorithm even 
as presently implemented , for two reasons. The first reason is that the ability to 
control zero-crossing times is not an absolute requirement for a pulse-shaping algo¬ 
rithm. The second reason is that the acceptable level of error for the VXIbus control 
strategy has not yet been defined. Of course, improvements in these two areas will 
make the algorithm even more useful. With respect to reducing noise, if the SNRs 
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of the TDW and of the RF pulse can be improved by 10 dB and 5 dB, respectively 
(’42), the algorithm will keep the pulse peaks in tolerance continuously. 

Third, power supply droop at dual-rated stations introduces only transient 
errors into the algorithm’s convergence. This causes the controller no significant 
problem. 

Fourth, a near-optimum TDW for the transmitter/antenna system can be 
obtained successfully off-line using the dummy load. In this way, when switching to 
the antenna, the newly designated operate transmitter comes on-line in tolerance or 
nearly in tolerance. 

B. RECOMMENDATIONS FOR FURTHER STUDY 

Further study is worthwhile in at least five areas. First, a reliable method to 
improve SNR for different bit resolutions will significantly increase the robustness 
and effectiveness of the steepest descent algorithm. Simple averaging, lowpass or 
bandpass filtering, and adaptive equalization are three possibilities. Second, other 
control algorithms may perform better than the steepest descent algorithm, partic¬ 
ularly if they are more robust and can control zero-crossing times automatically. 
Incorporating adaptive algorithms such as the recursive least squares method or the 
Kalman filter may work well. Third, a more effective strategy for controlling an 
entire PCI can possibly be found. For example, a better order in which to control 
the pulses might be pulse 1 (GRI A), pulse 1 (GRI B), pulse 2 (GRI A), pulse 2 
(GRI B), etc. Fourth, defining an acceptable level of error for the control process 
will be helpful. Keeping the pulse peaks in tolerance as defined by the seven tests 
in the signal specification for 100 percent of the control iterations may be neither 
practical nor desirable and may even be impossible. Perhaps an update to the sig¬ 
nal specification may become appropriate. Finally, writing a MATLAB function to 
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compute ECD for each iteration will provide statistical information on the effects of 
white noise and quantization error on ECD. 
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APPENDIX A 


COMPUTATIONAL METHOD OF ESTIMATING ECD 

(USCG Academy, New London, Connecticut) 

The General Problem 

The Loran-C antenna base current waveform can be expressed as 

x(t) = 0 ; t < t , 


and 


x (0 = A j ^-^ exp 

= Ar{t) , 


1 - 


65 J 


sinu 0 t 


where 

t is time in microseconds 

r is time origin of envelope (ECD) in microseconds 

A is pulse envelope peak in amperes 

u 0 is angular carrier frequency, 0.2 tt rad/^ sec 


The process of adjusting the TDW to establish an ECD and maintain some 
desired shape of the output pulse by visual comparison with a reference envelope 
(i.e., “pulse building”) can be thought of as a curve fitting process. 

The algorithm that is described accomplishes a MMSE fit that minimizes the 
squared difference between a set of eight half cycle amplitudes and some reference 
envelope of amplitude A and ECD r. The process of visually matching these two 
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data sets when expressed mathematically becomes a cost function, J. This squared 
error then becomes 


j = Ewo ~ Ar (*)l 2 . 

t =i 

where r(i) is the model which is a function of ECD. When J is minimized, this 
constitutes a MMSE fit. 


Minimization of the Cost Function, J 

In order to minimize J we will use partial derivatives. For well behaved Loran- 
C pulses and quadratic cost function, there is only one global minimum, and no 
maxima. Therefore we will set 


and 


Therefore, 


which implies 


and 


dJ_ 

dr 

8J 

dA 


= 0, 


= 0 . 


8 

I 

i=l 


J = EKO - Ar (0] 2 . 


= E 2 I 5 (0 - MOJl-K*)] = o , 


= £ 2 K*)-^r(i)] 


i=l 


-A 


r(i) 


= 0 . 


The solution for A is straightforward, yielding 


X>(*)r(0 


A = 


«=i 


E- 2 (0 

i=i 


The solution for r is a bit tougher! 
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Quadratic Approach 


J = £[«(*) - MO] 2 

1=1 

is called a quadratic cost function since for linear differences of [s(i) — Ar(i)], J is a 
second order polynomial. Although [s(i) — Ar(i)] is not a linear difference function 
of ECD, it becomes approximately linear in the region of minimum J for small 
differences of ECD. This says 

J - Jo- K(t - t 0 ) 2 , 

where 

Jo = minimum cost , 

and 

t 0 = the associated ECD at that J 0 

Now let’s choose three points for this function, separated by a common dis¬ 
tance, 8. This says 


(Jx-Jo) = K(t n - 8 - r 0 ) 2 ; t = t n - 8 

(Ji-Jq) = K{t n -t o ) 2 ; t = t n 

{J3-J0) = I<(t n + <5 - t 0 ) 2 ; t = t n + 8 
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However, the tq above does not provide an exact solution to the minimum of 
J. We’ll need an iterative algorithm. This algorithm can be stated as follows: 

Sk(J 3 ~ Ji) 

Tk+l " T * 2(Ji-2 J 2 + J 3 ) 

a) Let initial ECD = 0, Si = 1, compute Ji, J 2 , J 3 , T i 

b) Let 62 = 0.1, compute J \, J?, J 3 , t 3 

c) Let S 3 = 0.01, compute J u J 2) J 3 , U 

d) Let 63 = 0.001, compute J x , J 2 , J 3 , r 5 

t s represents the best estimate in the MMSE sense for the ECD. 
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APPENDIX B 

DATA COLLECTION AND FORMATTING 

The enclosed MATLAB programs describe how data vectors were collected and 
formatted for this project. These programs format the original ASCII data vectors 
and store them as vectors in eight sets in binary MAT-file form. This allows the 
data to be loaded quickly and easily. 
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% SAV_DAT: Corrects and formats Loran-C ASCII data files for use in 
% MATLAB & saves them in binary form in sets of manageable size. 

% In each set are several pairs of input and output vectors, labeled 

% xP and yP respectively, where P is the pair number, unique throughout 

% all the sets. Each set is a single MAT file named datasets.mat, 

% where S is the data set number. The variables may be loaded one 

% set at a time or all at once UBing LORLQAD. 

% Calls: MAX_VOLTS 

% Data pairs available: 

% 31 July 91: 5 pairs for the 42 xmtr with antenna simulator 

% 06 Sept 91: 3 pairs for the 42 xmtr with antenna simulator 

% 3 pairs for the 42 xmtr with dummy load 

% 28 Peb 92: 20 pairs for 42 xmtr w/ simulator 

% 20 pairs for 42 xmtr w/ dummy load 

% (Note: this data set is subdivided as follows: 

% 22 May 92: 16 pairs: an entire GRI, compensated and uncompensated 

% for droop (42 antenna simulator) 

% 30 Jun 92: 19 pairs (20, but pair 80 is actually 2 inputs) 

% for the 1 44A xmtr, with dummy load, antenna sim. 

% & 625' monopole antenna (pairs 81-83 only are 

% on antenna) 

% This program decimates each data vector to a desired sampling 
% frequency, using a lowpass filter to prevent aliasing. The zero 

% level of each data vector, estimated by taking the mean of the 

% last 596 sample points, is subtracted so that any DC measurement 
% bias is removed. Bach data vector is normalized to 1 and then 

% then scaled up to the correct voltage. Although the Loran-C output 

% current is customarily measured, in this engineering model the output 
% voltage is measured, using an infinite input impedance at the 
% oscilloscope. This is not the same as the voltage read by the 

% Electronic Pulse Analyzer (SPA) . In cases of 6 Sept 91 and following, 

% the transmitter cathode current (TKI) was generally held at about 1.0 
% amp. All ’42 rf vectors were measured from the J6 jack. The '42 
% rf vectors were measured from the J26 jack. For the ! 44A dummy 

% load data a 20 db attenuator was used (model 42, ser # 173-56) 

% since the J26 jack is uncalibrated and unloaded. If implemented 
% on the 1 44A this should really be buffered better. 

% Dean C. Bruckner, 11/22/91 Rev. 9/4/92 

clc; 

dispC Caution: running this program will clear the workspace.'); 
dispC Press return to continue or ctrl-C to abort.') ;disp('') ; 

pause 

clear 

Nl» [1 2 4 8] ; 

Nl_indx>menu (' Select Decimation factor ', ' 1',' 2',' 4',' 8') ; 

JUKI(Nl_indx);disp{['N * ',num2str(N),' selected’]) 
dispC Press <Bnter> to continue or ctrl-C to abort'), pause 
fs-10e6/N;len-409S/K; 


pt«[l 6 12 22 32 42 
5 11 21 31 41 51]; 


% Data pair numbers in each set 








[xx,nsets]-size(pt);clear xx; 
max volts 


% Load max volts and dc bias 
% from o-scope measurements 


for s-l:neats; 

S-num2str(s) 
x_str«' ';y_str«' '; 

n-1; % Loop index 

for p-pt(l,s):pt(2,s) 

P-num2str(p) 

aval(['load plot',P,'a.dat; load plot',P,'b.dat;']) 

eval(['x«plot',P,'a(l:4096);']) 

xsx' -mean(x(3501:4096)); 

aval(['y-plot*,P,'b(l:4096);']) 

y»y'-mean(y(3501:4096)); 

X-fft(x,8192);Y-fft(y,6192); * Apply ideal LP filter 

startf-8192/(2*N)+l;endf-8192-startf+l; 

X(startf :endf) * [] ;Y(startf :endf) - [] ; 
x-real(ifft(X));y»real(ifft(Y)); 
x(startf :2* (startf-1) )*[]; 
y(startf:2*(startf-1))«[]; 

x«x/max(x) * (mv_in(n,s)-dc_in(n,s)); 
y»y/max(y) * (mv_out(n,s)-dc_out(n,s)); 

plot (x) ;pause (2) ,-plot (y) ;pause (2) ; 
aval (['x' ,P, ' «x;' ]) 
aval (['y',P, '=y;']) 

aval(['clear plot' ,P, 'a plot',P,'b']); 

x__str» [x_str, ' x', P] ;y_str» [y_str, ' y', P] ; 
n-n+1; 
end 

if s*«l; 

y4t*y4; % correct synch, error 

y4(1:48/N)*y4t(4096/N-48/N+1:4096/N); 
y4(48/N+l:4096/N)*y4t(1:4096/N-48/N); 

diap('y4 adjusted to correct synchronization error');disp(’’); 
and 

aval(['save dataset',S,'.mat ',x_str,y_str]) 
disp(['dataset',S,'.mat:']); 
disp([x_str]);disp([y_str]); 

and 

% Droop data set for '42 (set 7) & '44 data (set 8) 
pt* (52 68 
67 87]; 

[xx.nsets]-size(pt);clear xx; 
for s*l:nsets; 

num_pairs-pt(2,s)-pt(l,s)+l; 
if s.*l 

x_scale»4*.2*ones(l,num_pairs); % .2V/div (.8V 0-pk) 

y_scale*4* 5*ones(l,numjpairs); % 5V/div (20V 0-pk) 

alseif s--2 % Here use these to invert also 

x scale-4*[.5*ones(1,8) 1122 127/2 -.5 .5 -1 .5*ones(1,4)]; 
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127/2 2-22 -2*ones(1,4)]; 


y_scale»4*[-2*ones(1,12) 
end~ 

S*num2str(s+6) 
x_str«' ';y_str«' '; 

n»l; ~ % Loop index 

for p*pt(l,e):pt(2 ,b) 

P«num2str(p) 

eval(['load plot',P,'a.dat; load plot',P,'b.dat;']) 

aval(['x-plot',P,'a(1:4096);']) 

x«x'-mean(x(3501:4096)); 

eval(['y-plot',P,'b(l:4096);']) 

y«y'-mean(y(3501:4096)); 

mx« m ax(x);my«max(y); % Save to restore orig. scaling. 

X«fft(x,8192);Y«£ft(y,8192); % Apply ideal LP filter 

etartf-8192/(2*N)+1,-endf-8l92-startf+1; 

X(startf :endf)« [] ;Y (startf :endf)»[]; 
x«real (ifft (X)) ;y»real (ifft (Y)) ; 
x(startf:2*(startf-1))«[]; 
y(startf:2*(startf-1))«[1; 

xax/max(x)*mx/l27*x_scale(n); % Scaled differently than 

y»y/max(y) *my/127*y~scale(n) ,* % sets 1-6. 8-bit resolution 

% used in Lecroix o-scope 

if gasl 

if any(find(p«»[55 57 67 ])) 

y-y; 

else 

y» -y; % Correct invertion introduced 

% at J6 jack (except for 

% incorrectly sampled vectors, 

% where correction would make 

% phase code incorrect again). 

end 

if any(find(p>s65)) % Correct incorrectly sampled 

x« -x; % input data vector (i.e., 

end % pos phase GRI captured 

end % instead of negative) 

plot(x);text(.15,.16,P,'sc'),pause; 
plot(y);text(.15,.16,P,'sc').pause; 
eval(['x',P,'*x;']) 
eval([’y',P,'*y;’]) 

eval(['clear plot',P,'a plot',P,'b']); 

x^str-[x_str,' x',P];y_str«[y_str,' y',P]; 
n-n+1; " 

end 

eval([’save dataset',S,'.mat ',x_str,y_str]) 
disp(['datasetS,'.mat:']); 
disp([x_str]);disp([y_strj); 

end 
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% MAX_VOLTS: Measurements from oscilloscope plots for AN/FPN-42 
% Loren-C transmitter. Bach value if the max positive value 
% in volts. Called from SAV_DAT. Data sets 7 and 8 were 

% formatted more directly and easily by just using the o-scope scale 

% instead of reading the plots. But since these were done already 
% they stayed as is. 

fflv_in«zeros(10,nsets); 
mv_out-zeros(10,nsets); 
dc_in-zeros(10,nsets); 
dc_out-zeros(10,nsets); 

mv_in(1:5,1)« [.72 .74 .28 .28 .48]'; % (from o-scope plots) 

mv out(1:5,1)- [17.2 16.2 3.3 1.62 17.2]'; 

mv”in(l:6,2)- [.82 .8 .93 .9 .52 .51]'; 

mv_out(1:6,2)- [24.5 8.2 27 10 3.0 1.2]'; 

mv in(:,3)- 2*[2 .85 .6 .55 .5 .9 3.3 2.3 1.8 1.35]'; 

nrv — out (:, 3)» -2*10*[1.S5 1.45 1.15 1 .83 .5 .26 .22 .2 .27]'; 

mv”in(:,4)« 2*[4 4 2.75 .86 .75 .63 .46 .3 .32 .4]'; 

mv”out(:,4)» -2*10* [1.4 1.37 1.5 1.37 1.29 1.2 .68 -1.3 -1.39 -1.65]'; 

mv”in(:,5)« 2*[2.4 1.05 .69 .61 .49 1.9 3.75 2.6 1.1 1.09]'; 

mv”out(:,5)- -2*10*[.49 .47 .39 .265 .22 .46 .14 .12 .135 .141]'; 

mv_in(:,6)» 2*[3.65 3.95 3.79 1.25 .9 .8 .79 .305 .36 .54]’; 

mv”out(:,6)« -2*10*[.2 .46 .47 .46 .43 .38 .39 -.4 -.44 -.53]'; 

♦ Corrections: 

% - Outputs accidentally inverted 

% 2 Correcting for differences in o-scope input impedance 

% settings (50 ohm evenly divided voltage with the 

% load and so its amplitude was only half of those 

% taken with infinite input impedance. Thus they need 

% to be doubled. 

% 10 A 10X probe was apparently left out of these meas. 

dc_in(l,l)» .01; % DC voltage bias in measurements, 

dc out(1:5,1)- [.1 0 .05 .02 0]'; % from plots 

dc”out(1:6,2)- [.5 .2 .2 .2 .05 0]'; 

dc_in(:,3)- .01*[3 zeros(l,8) -1]'; 

dc”out(:,3)« .01* [3 zeros(1,5) .5 0 0 0] ' ; 

dc_in<:,4)- .01*[zeros(1,5) .5011 1]'; 

dc_out(:,4)» .01* [zeros d,7) 1 2 2]'; 

dc in(:,5)— ,01*[3 211155552]'; 

dc“out(:,5)« .01* (0 0 1 1 1.5 1 .5 .5 .5 .5]'; 

dc~in(:,6)« .01*[10 55333311 1]'; 

dc OUt(:,6)- .01* [.5111111211]'; 


% Max volts of input vectors 
% Max volts of output vectors 
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APPENDIX C 


FITTED CURVES FOR CATENATED MODELS 


This appendix contains fitted curves of magnitude and phase for the catenated 
models, as follows. For values of E n outside the endpoints of the curve, the function 
simply assigns the value of the curve at the nearest endpoint. In the titles “A” 
indicates the antenna and “D” the dummy load. 


Transmitter 

Description 

Page 

AN/FPN-42* 

19-point curves for combined model, antenna 
simulator (4 poles, 3 zeros) 

148 

AN/FPN-42* 

13-point curves for combined model, dummy 
load (2 poles, 1 zero) 

152 

AN/FPN-44A* 

15-point curves for combined model, antenna 
simulator (6 poles, 5 zeros) 

154 

AN/FPN-44A* 

4-point curves for combined model, dummy 
load (4 poles, 3 zeros) 

160 

AN/FPN-42 

5-point curves for first catenated model, an¬ 
tenna simulator only (4 poles, 3 zeros) 

164 


* 4-part set of curves used in computer simulation program for N = 1. 
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CPLX POLE PAIR 2, ’42 (A, N = l) 


Magnitude of Solaetad Parameter 
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REAL ZERO, ’42 (A, N=l) 


Mognitud* of S«l«et«d Poram«t«r 
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REAL ZERO, ’42 (D, N = l) 


Magnitude of Selected Parameter 
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CPLX POLE PAIR 2, ’44A (A, N = l) 



Pha*« of S«l«et«d Porom«t«r 
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CPLX POLE PAIR 3, ’44A (A, N = l) 



Phase of Selected Parameter 
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CPLX ZERO PAIR 1, ’44A (A, N = l) 
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Energy of normalized input vector, wott-eec X ie-6 
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CPLX ZERO PAIR 2, ’44A (A, N=l) 




Phoee of Selected Parameter 



Energy of normalized Input vector, wctt-»ec X le-6 


























REAL ZERO, ’44A (A, N=l) 
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CPLX ZERO PAIR 1, ’44A (D, N=l) 
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CPLX POLE PAIR 1, ’42 (A, N=l, 5 PAIRS) 
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CPLX POLE PAIR 2, ’42 (A, N = l, 5 PAIRS) 


Mo$nitud« ef S«l«et«d PoromoUr 



Pho$# of S*l«Ct«d Porom«t«r 
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CPLX ZERO PAIR, ’42 (A, N=l, 5 PAIRS) 


Mognitud* of s«l«ct«d Porom#t*r 



PhoM of S«l*ct«d Porom«t*r 



Entrgy of normollxod input vector, woU-*#c X 1*-6 























REAL ZERO, ’42 (A, N=l, 5 PAIRS) 


Magnitude of Selected Parameter 



Energy of normalised input vector, wott-eec X 1e-6 
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APPENDIX D 

SELECTED MATLAB PROGRAM LISTING 


The following MATLAB M-files are part of the simulation program. The main file 
is SIMZ, which calls more than fifty functions in the course of a simulation run. 
Space does not permit a full listing here. The last function listed, MODEL_YW, 
was used for modeling and is not directly part of the simulation program. 
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% SIMZ: Simulates controlling the AN/FPN-42 and 44A Loran-C 
% transmitters in the z domain. 

% Calls SETUP, AR2, XMTR, INTERP, ERRORS, PULS_INI, MAINMENU, 

% FIND_PSV, ENVEL, DISPZ 

% and the function defined in 'fcall' (see SETUP). 

% Also calls file named in row string matrix ALG_SWP (See SETUP) 

% Deem C. Bruckner, 4/22/92, rev. 9/12/92 

% ***************** initialize program *********************** 

clear 
setup 

first_run='y'/ 
done_simz='n'; 
mainmenu 

while done_simz=='n' 

If **************** Enter control loop ************************* 
for p=l:length(control)/ 

if p**l,-num_iter*num_iterl/else num_iter=num_iter2;end 
% If drift will be used, start with nonzero random drift vector 
% if first run. 

if drift_iter~*0 & first_run=='y'/ 

[drift,drift_stor]*ar2(drift_stor,step); 
end 

PSTR=['Pulse ',num2str(control(p))]/ 

PC=phasecode(control(p))/ 

tdw*tdw_pci(:,p)/ % Select one input pulse 

puls_ini % Initialize pulse convergence 

% plot and misc error matrices 

% for this loop only. 

done jpulse*'n';m=l;md=l;ms=l/ 
tOsdock; 

while done_pulse=='n' 

if m>l/times(m)setime(clock,tO)/end 
tOsdock/ 

M=num2str(m)/Ml=num2str(m-1)/ 

% ************ Find resti for GRI of subject pulse *************** 

restl*rests(rest_indx)/ % Rest time for pulse 1 of GRI 

blank*blanks(rest_indx) ,- % l=rate blanked in this iter. 

% 0=rate not blanked. 

rest_indx*rest_indx+T_dso/ratel/ % Find rest_indx for next time, 

while rest_indx>length(rests)-1 % Wrap index around 

rest_indx*rest_indx-length(rests)+2/ 
end/ 

if m<5 & err_disp=-l/plot(tdw)/pause/end 
if blank-*!/ 


V Declare vars & intialize 

% First menu (see end of loop 
% for other occurrence) 












if err_disp»«l; % Skip blanked GRI 

disp([ 1 Rate 1 blanked between m = 1 ,Ml, 1 and m = ',M]) 
end 

blanksav-[blanksav m]; 
else 

% ********** Get ps_volt for subject pulse in GRI ************** 

ps_volt-find_psv(restl)*psv_sim(control(p)); % Find rest time 

% for beginning of new GRI 
V & estimate ps_volt for 

% p-th pulse 

% ********** Generate & capture rf ***************************** 

if step_iter~*0 & ms>step_iter & p=«l; % If time has arrived, 

step-1 ,-m_step-[m_step m]; % switch transmitters & record 

ms*l ; 

if err_disp«=2; % Inform user 

text7m,4.3e-4,'x');text(m,4.37e-4,' A '); 

else 

disp(['Switched transmitters at iteration ',M]) 
end 
else 
step-0; 
ms=ms+l; 
end 

if drift_iter *"« 0 & (md>=drift_iter|step==l) 

[drift,drift_stor]-ar2(drift_stor,step); 
md*l ; 
else 

md»md+l; 
end 

rf*xmtr(tdw,PC,drift,ps_volt); 

if m<5 & err dispssi ,-plot (rf) ,pause;end 

if m»l | skTp_flag==l;y*envel (rf ,tdw, PC) ;end 

% ********** Calculate, display & record errors **************** 


if m««l 

err_sav(m,:)^errors(y,yO,PC,m,err_disp); 

else 

err_sav(m,:)-errors(y,yO,PC,m,err_disp,err sav(m-l,:)); 

end 

if any(err_sav(m, 2:4)>[.01 .03 .1])==1 

OOT(m)=1; % Record out of toler. iterations 

end 

if err_disp*=2; % Flash if out of tolerance 

semilogy(num_iter/9,1.3e-4,'w*'),hold on 
if OOT(m)*=1 & round(m/2)--m/2 

semi logy(num_iter/9,1.3e-4,'i**).hold on 
end % Note: 'round' used to ident. 

end % every other iteration. 

ps_sav(m,:)-[restl ps_volt]; 
drift_sav(:,m)-drift;” 
y_savT:,m)-y; 
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%snr_tdw(m)«snr(tdw ); 
%snr_r£ (m) >snr (rf) ; 


% record SNRs for averaging 
% Note: change DISPZ also 


I; ********** produc6 nsv tdw *********************************** 
eval(f_call) 

% ********** Pill AFG buffer with new inputs ******************* 

for pp*p:length(control) % (pth & following pulses) 

% If in GRI A, boost each following tdw in GRI A so when it is 
% controlled, convergence will be faster. Undo the phasecode 
% of each pulse & then reapply it as appropriate. Variable 

% "boost" is set in ****_INI, not in XMTR_CFG (since knowledge 

% of needed boosts should be experimentally obtained). 

% To boost the tdws in GRI B when controlling a pulse in GRI A, 

% scale back to pulse 1 in GRI A and then skip to GRI B, scaling 

% iip from the first pulse in GRI B. 

% If in GRI B, do the same. 

if control(pp)<=len_p/2; % both in GRI A 

tdw_pci (: ,pp) -boost'*' (control (pp) - control (p)) *. . . 

(tdw*(-1)“PC)*(-1)“phasecode(control(pp)); 
elseif control(p)<-lenjp/2 % in different GRIs 

tdw_pci(:,pp)»boost“{1-control(p))*_ 

boost“(control(pp)-lenjp/2-1)*... 

(tdw*(-1)“PC)*(-1)“phasecode(control(pp)); 
else % both in GRI B 

tdw_pci(:,pp)=boost“(control(pp)-control(p))*... 

(tdw*(-1)“PC)*(-1)“phasecode(control(pp)); 

end 

end 

% ********** swap XMTR loads if error below threshold ********** 

if xmtr_load== 1 Dummy Load 1 & all(err_sav(m,2:4)<[.006 .015 .05])==1 
m_swp=m; 

xmtr_load = 'Antenna '; % Change to antenna 

eval(alg_swp(alg,:)) % Reset part of algorithm 

if err_disp=*2; % Inform user 

text(m,4.3e-4,'a');text(m,4.37e-4,'“'); 
else 

disp(['Switched to antenna after iteration ’,M]) 
end 

skip_flag=l; 
end 

m*m+l ; 

end % end 1 iteration 

if m>«num_iter+l;done_pulse='y';end 

end % end 1 pulse 

% ************** save rf and display results for pth pulse ********* 

rf jpci(:,p)=rf; 
if err_disp**2 
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text(.47,.2,'Press <Enter> to continue','sc'),pause 
hold off 
else 

dispC Press <Enter> to continue') ,pause 
end 

if any(control(p)==analyze)«*1;dispz;end 
end % end 1 simulation run (all puls) 

first_run='n' ,- 

mainmenu % get new parameters for more runs 

end % end all runs & leave SIMZ 


0 


% SETUP: Sets up workspace for Loran-C simulation file SIMZ. 

% Declares global variables and assigns values to them. Sets up 
% workspace format and declares default settings such as the type 

% of random number distribution, etc. When user selects a control 

% method a script file is used to initialize only those control 
% variables for the selected algorithm, including the control 

% function call in a text string. Note that when a default setting 

% is changed here, it should also be changed in MAINMENU. If 
% an algorithm is added, review MAINMENU carefully to ensure 
%• parallel changes are made properly. 

% Calls PGEN, FLIPLR, ENVEL, and files named in string matrix ALG_INI 
% ... COEFFSAV.MAT (created by M0D1FMT1) , H_SCRPT, EOUTJEIN.MAT, 

% Declares all global variables used in this simulation. 

% For complete list of variables see INDEX. 

% Dean C. Bruckner, 4/22/92, rev. 9/2/92 

format compact 
format short 
rand (' normal') 

clc;disp('Initializing ...') 

% **************** INITIALIZE TRANSMITTER MODEL ********************** 

% **************** Declare global variables 

global phasecode ratel rate2 sta_id 

global coeff bound root_indx ps_id cp cn cc cr 

global drift_sdev drift_ref ar2_var ar2_a ztoler ztimes 

global eout_ein ps_tau ps_prev ps_imp ps_indx N fs fs_adjust 

global xmtr_id xmtr_load bits imbalance sig_noise yOD yOA tdwjnoise 


% **************** s e t U p dual rate parameters 

load restvars % Load variables for dual-rate 

% simulation (rests, sta_id, 

% ratel, rate2) 

T_dso*4; % Interval at which rf pulses 

% are sampled by the digital 
% storing oscilloscope (secs) 


rs<*rand( 'dist') ,-rand( 'uniform') % 

rest_indx*round(rand*.9*length(rests)); % Start rest times index randomly. 
rand(rs) 
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t 


Set up Pulse measurement parameters 


stimess(-25:5:30)'*le3; % Zero crossing times (ns) 

ztoler=[ 2000 1500 1000 500 250 0 100 100 100 100 100 100 

1000 100 75 50 50 0 50 50 50 50 50 50] ' ; 

% These correspond to the category 2 & 1 xmtrs in the Loran signal spec. 
% The *42 is a cat. 2 xmtr; the '44A is catl (Note that xmtr_id is 
t opposite from this: xmtr_id*l is for the *42; xmtr_id=2 is the ’44A) 


% **************** Configure transmitter 

xmtr_id*l; % 1 * 42, 2 = 44A 


xmtr_load='Dummy Load'; % Starting load of xmtr, 

% 'Antenna ' or 

% 'Dummy Load' (note: string 

% lengths must be equal!!) 

if sta id(1)==' S'; % Phasecode: 0=pos, l=neg 

phasecode*[0000 0110 0101 0 0 11']; 

else 

phasecode^[0011 0101 0 0110 0000 1' ] ; 

end 

len_p«length(phasecode); 


imbalance*0; % [0-100]; pet. by which the ampl. 

% of neg. phase coded pulses is 
% decreased. 

% **************** Load curves of xmtr model 

N*l; % Default decimation factor 

xmtr_cfg % Configure transmitter 

% **************** s et up parameters for simulation 

skip_flag*0; % Used in switching automatically 

% from dummy load from antenna 
control=l; % Pulses to control 

analyzesl; V Pulses to analyze (after 

% convergence) 

tau*0; % Assigned local BCD 

ETA*0; % Barly timing adjuBt in 

% microsecs (used in PGEN) 

bitSoS; % Function generator resolution 

% ('O' selects best floating 

% point resol. of the computer) 

num_iterl*100; % Iterations of 1st selected pulse 

num_iter2=20; V Iterations of following pulses 

step_iter*0; % Interval between xmtr steps 

% (switches) 

err_disp=2; % Method of displaying errors 

% during convergence 
% 0=none,l=text,2=plot 


% ****************** INITIALIZE CONTROL ALGORITHM ******************** 

alg*l; % Default algorithm to start 

% The following 5 string matrices are used to handle the current 
% algorithm without listing the names of the algorithms or their 

% associated files anywhere else in the program than here. The first 

% string matrix holds the names of the algorithms; the following 4 
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matrices hold names of files associated with each algorithm which 
are called at different points in the program. Lines 2 and following 
of each matrix should be the same length as line l and should follow 
suit. 


Example: 

alg_name=['Steepest Descent '; 

'Neural Network 1 ; 

1 Recursive LS ']; 

The Steepest Descent Algorithm uses 6 modular M-files: 

DESC The function 

DBSC_INI Initializes the algorithm 

DESCMENU Lets the user change parameters easily 

DBSCHEAD Menu header for DESCMENU 

DBSC_SWP Resets part of the algorithm when xmtr load is swapped 

DESCDISP Displays alg params for most recent run in DISPZ 

Other algorithms should use the same file structure. Details of the 
minimum requirements for each of the above files are listed in the 
text of each file. 


alg_name=['Steepest Descent 
alg_menu=[ 1 descmenu' ] ; 
alg_swp *['descswp 1 ]; 
alg_ini *['desc_ini']; 
alg_disp=['descdisp 1 ] ; 


' ] ; % Algorithm names 

t M-file names 

^ it n 

% n n 

% n n 


eval(alg_ini(alg,:)) 


% Initialize algorithm 



function rf=xmtr(tdw,PC,drift,ps_volt) 


Punction rf=XMTR(tdw,PC,drift,ps_volt): Simulates the AN/FPN-42 
or AN/FPN-44A Loran C transmitter. To account for nonlinearities, 
the poles & zeros of the xmtr's transfer flinction are modeled 
as a function of the normalized power supply voltage and the 
energy of the normalized input waveform. 

Uses global variables: bits, imbalance, sig_noise, xmtr_load 
Calls: ENERGY, FIND_AB, FIND_PSV 


Local variables: 
A 

a,b 

cr 

drift 

energy_in 

energy_no 

h 

lc 

ps_volt 

restl 

tdw 

xmtr_load 

rf 


Amplitude of output vector (volts) 

Denominator & numerator polynomials of model 
Number of rows in 'coeff' 

Parameter vector modeling xmtr drift 
Energy of input vector (watt -sec); R=l. 

Energy of normalized input vector (watt-sec); R=l. 

Transmitter inpulse response sequence 

Slope & intercept of energy in/out of xmtr 

Estimated normalized power supply voltage: (0,1] 

Power supply recovery time since last pulse 

Transmitter drive waveform 

String: Defines load connected to xmtr: 

'A'=antenna, 'D'=dummy load 
(Radio freq) output pulse 


Dean C. Bruckner, 7/17/92, rev. 9/7/92 


Obtain xmtr transfer function ********************* 
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if nargin<4;ps_volt=l;end V 

if nargin<3;drift*zeros (cr, 1) ,-end % 

if nargin<2;PC-0;end % 

len_tdw«length(tdw); % 

[energy in,energy_no]-energy(tdw); \ 

% 


[a,b]-find_abl(energy_no,ps_volt,drift); 
% **************** Produce output ****** 
rf-filter(b',a’,tdw); 

% **************** Adj us t output ******* 


if xmtr_load«' Dummy Load 1 ; % 

lc«eout_ein(2,:); % 

else % 

lc-eout_ein(1,:); 
end 

energy_out-lc(1)*energy_in + lc(2); 
rf*rf/max(rf); % 

A-sqrt(energy_out/((rf'*rf)*N*100e-9)); % 

% 

rf«A*rf*ps volt; % 

% 

% 

rf-rf-mean(rf); % 


Default: fully recovered 
Default: no drift 
Default: pos phase code 

Length of input vector 
Energy of input vector 
(regular & normalized) 

% Denom & num polynomials 


Apply input energy vs output 
energy to find output 
aaqplitude. 


Normalize output sequence & 
calculate final normalization 
power (power norm, to R-l) 
Assign estimated output energy 
to output sequence, including 
power supply droop. 

Remove DC bias 


% Note on tube xmtr imbalance: 

% Initially the xmtr imbalance (between the 2 banks of tubes which 
% anqplify the positive & negative parts of the pulse, respectively) 

% was modeled in detail, as shown below, by adding up to one percent 

% distortion to the positive samples. 

%rf(.ind(rf>0))-(1-.01*imbalance) * rf(find(rf>0)); % xmtr imbalance 

% Apparently this is an accurate representation of the distortion. 

% However, I could not find real documentation on the phase code 
% balance adjustment that described exactly how this was remedied, 

% just that the imbalance caused negatively phase coded pulses to 

% be smaller in amplitude (both pos & negative parts equally), as 

% described in LCDR Taggart's EERP & VXIbus reports. According to 
% him, the phase code balance adjustment simply increases the amplitudes 
% of the TDWs for the negatively phase coded pulses, not adding any 
% DC bias level, etc. Therefore, the imbalance is now modeled 

% as a percentage decrease in the amplitude of rf. This will 

% be compensated for automatically just like droop, since each pulse 
% of the PCI is controlled independently. 

if PC--1; % For negatively phase coded 

rf*rf*(1-.01*imbalance); % pulses, decrease amplitude 

end % by a percentage. 

if sig_noise~*0 % Misc white noise in output 

rs«rand('dist');rand('normal 1 ) % (std dev expressed as 

% percentage of peak anqpl.) 
rf-rf + .01*sig_noise*max(rf)*rand(length(rf),1); 
rand(rs) ; 
end 


if bits ”-0; 


% Amplitude resolution in DSO 
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max_rf*max(rf); 

rf *round(rf/ma.x_rf * (2 A bits)/2) / (2 A bits) *2*max_rf ; 
end 



function tdw*pgen(x,PC,ETA); 


Function tdw*PGEN(x,PC,ETA): Simulates analog pulBe generator, 
given a 16-element vector of peak voltage values (pos & neg 
or all pos). Pulse is triggered 10 us after beginning of 
data vector. Due to the problems of dealing with fractional 
values of samples per period, tdw is formed at 5.0 MHz and 
decimated down to the desired sampling frequency. 

Calls global variables: bits, xmtr_id, fs 


Variables: 
ETA 


fsjpg 

len_pg 

PC 

s_eta 

spp 

tdw 

x 


Early timing adjust, in microsecs. Changes phase 
of tdw within window (pulse still begins exactly 
at the trigger--ref. discussion w/ LCDR G. Kmiecik 
on 5/22/92) . The effect of the ETA shows up in 
the Envelope-to-Cycle-Difference (ECD) of the output. 
Sampling frequency used to build tdw 
Number of sanples in tdw at fs_pg 
TDW phase code (0=pos, l=neg) 

Number of sanples in ETA 

Sanples per period (period = 10 us) 

Transmitter drive waveform 

16-element input vector of pk voltage values 
(pos & neg or all pos) 


Dean C. Bruckner, 2/21/92, rev. 9/7/92 


******************* verify inputs ******************************** 


if N>*2 

fs_pgj=5e6;lenjpg=2048; % 5 MHz has a whole number of 

% sanples in each half-cycle. 

else 

fs_pg=fS;len_pg=4096; 
end 

if nargin<3;ETA=0;end 

if nargin<2;PC*0; % Check phase code 

elseif PC~*0 & PC“*1;error('Phase code must equal 0 or 1'); 
end 

if nargin<l,-x*ones (16,1) ; % Default: all 1/2 cycles equal 

elseif size(x)**[1,16];x=x'; % Reorient if necessary 

end 

if size(x)“*[16,1];error('Size of x incorrect');end 
x«abs(x); 

% ************ Generate PGEN input vector ************************** 

s_eta*round(ETA*le-€*fsjpg); 
spp*1Oe-6*fs_pg; 
k»(l:8*spp) 1 ; 
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m-ones(spp/2,1)*x'; 

m* reshape(m,round(8*spp), 1) ; 


% Extend vector x to each "bin 
% which is mult, by tdw 


tdw(1:spp, 1) =zeros (spp, 1) ; % Trigger pulse at 10 us point 

templ»sin(2*pi*k/spp - pi*PC). *m; V Generate sine pulse 

teiqpl (length (tempi)-s_eta+l: length (tempi) ) = [] ; % Apply ETA 

tenqp2*-flipud(tempi (l:s_eta)) ; 

tdw (spp+1:9*spp) ■= [temp2; tempi] ; % Recombine 

if xmtr_id==2 

nn=l:len_pg-9*spp; % Tail drive circuit 

tdw(9*spp+l :len__pg> =x(16) *. . . 

(exp((-(le€*nn/fs_pg)/22)}. *sin(.2*pi*le6*nn/fs_pg-pi*PC))'; 

else 

tdw(9*spp+lzlenjpg)=zeros(len_pg-9*spp,1); 
end 

if N>2 

tdwatdw(1:round(fs_pg/fs):len_pg); % Decimate (ignore aliasing) 

end 

if bits "=0; % Amplitude resolution in AFG 

max_tdw*max(tdw); 

tdw»round(tdw/m«uc_tdw* (2*bits)/2) / (2'‘bits) *2*msoc_tdw; 
end 

if sig_noise~=0 % Misc white noise in input 

rsz=rand( 1 dist') ;rand( 'normal') V (std dev expressed as a 

% percentage of peak amplitude) 
tdw»tdw + .01*tdw_noise*sig_noise*max(tdw)*rand(length(tdw),1); 
rand(rs); 
end 



function ps_volt=find_psv(restl) 

% Function ps_volt*FIND_PSV(restl): Given the time the xmtr power 
% supply has had to recover from the last pulse of the preceding 
% GRI, this function estimates the new power supply voltage 
% (norialized and in the range (0,1]) for pulse 1 in the new GRI. 

% Uses global vars: ps_tau ps_prev ps_imp 
% Calls: 

% Variables: 

% ps_volt Estimated normalized power supply voltage: (0,1] 

% restl Power supply recovery time since last pulse 

% Dean C. Bruckner, 7/20/92, rev. 9/7/92 

if restl<.001;error('restl must be >= .001 sec.'),-end 
t«-ps_tau*log(1-(ps_prev-ps_iiqp)); % Point on the curve where 

% recovery starts at end of 
% last GRI (note: "log" is 
% the natural logarithm) 

ps_volt»l-exp(-(t+restl)/ps_tau); % ps_volt after resting "restl" 

% seconds 


178 









function [a,b,h_yw] =model_yw(h,p,q) 


% Function [a,b,h_yw]=MODEL_YW{h,p,q): Solves Yule-Walker equations to 
% find pole-zero model of Loran-C data vector. 

% Dean C. Bruckner, 4/7/92. Adapted from algorithms written by 
% Tom Johnson of the Naval Postgraduate School. 

% Ref: C. W. Therrien, Discrete Signal Proc. & Statistical Signal 

% Processing, Prentice-Hall, 1992. 


len_h*length(h) ,- 
Rx»lorcorr (h) ,- 

R»toeplitz(Rx(1:len_h),Rx(1:p+l)); 
Re»R(q+2:len_h,:); 

[re, ce] *size (Re) ; 

a= [1; -Re (:, 2 : ce) \Re (:, 1) ] ; 


% Autocorrelation vector by ffts 
% ACF mtx 

% Extended corr mtx (Ther. 9-168) 


% using extended Y-W method 

% Note: Matrix pinv calculated 

% by gaussian elimination. 


ha=filter(l,a,[1;zeros(len_h-l,1)]),- 
Ha=toeplitz(ha,[ha(l),zeros(l,q)]); 
b=inv(Ha'*Ha)*Ha'*h; 


% Use Shank' s method to 
% find "b" parameters. 
% (Therrien ch. 9) 


h_yw=filter(b,a,[l;zeros(len_h-l,1)]); 


% ARMA model impulse resp. 


disp('Poles: '); 

[abs(roots(a)) angle(roots(a))] 
disp('Zeros: '); 

[abs(roots(b)) angle(roots(b))] 

disp('Press <Enter> to show plots');pause 

clg 

if q>0;rts_b=roots(b); rts_b=rts_b(find( abs(rts_b)<2 )); 

polar (angle (rts_b) ,abs(rts_b) , 'go') ;hold on,-end; 

polar(angle(roots(a)),abs(roots(a)),'rx'); 

grid;title('Pole/zero plot of the Yule-Walker estimate'); 

pause 

hold off 


clg;plot (1 .-length (h) ,h, ' - ', 1 .-length (h_yw) , h_yw, ' - - ') 

grid,-title ('Original & modeled sequence') 

xlabel (' Sample number ’) ,-ylabel (' Magnitude') 

text(.5,.15,['MSB = ',num2str(mse(h,h_yw)*le6),'e-6'],'sc') 

pause 
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